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SUMMARY 

The  long  and  intermediate  length  Josephson  tunnel  junction  oscil¬ 
lator  with  overlap  geometry  of  linear  and  circular  configuration,  is 
investigated  by  computational  solution  of  the  perturbed  sine-Gordon 
equation  model  and  by  experimental  measurements.  The  model  predicts  the 
experimental  results  very  well.  Line  oscillators  as  well  as  ring  oscil¬ 
lators  are  treated.  For  long  junctions  soliton  perturbation  methods  are 
developed  and  turn  out  to  be  efficient  prediction  tools,  also  providing 
physical  understanding  of  the  dynamics  of  the  oscillator.  For  inter¬ 
mediate  length  junctions  expansions  in  terms  of  linear  cavity  modes 
reduce  computational  costs. 

The  narrow  linewidth  of  the  electromagnetic  radiation  (typically 
1  kHz  of  a  line  at  10  GHz)  is  demonstrated  experimentally.  Corresponding 
computer  simulations  requiring  a  relative  accuracy  of  less  than  10-7  are 
performed  on  supercomputer  CRAY-1 -S.  The  broadening  of  linewidth  due  to 
external  microwave  radiation  and  internal  thermal  noise  is  determined. 


The  effect  of  constant  magnetic  fields,  applied  to  tune  the  radia¬ 
tion  frequency,  on  the  resonant  soliton  oscillations  is  investigated  by 
detailed  computations  of  the  power  spectra.  Hysteresis  and  chaotic  inter- 
mi  ttency  between  soliton  dynamic  states  are  found.  In  narrow  windows  of 
parameter  space  chaos  effects  cause  noise  rise. 

Poincarg  and  return  maps,  Painlevg  and  Melnikov  methods  are  applied 
to  indicate  and  predict  chaos  in  ordinary  differential  equations  of 
sine-Gordon  type  modelling  SQUIDs. 


KEY  WORDS 


Cavity  mode.  Chaos,  Fiske  steps.  Hysteresis,  Intermittency , 
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tunnel  junction.  Noise  rise.  Perturbation  theory,  Radiation  line- 
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In  the  following  sections  we  summarize  the  new  main  results  ob¬ 
tained  under  the  contract.  References  to  the  published  papers  in  the 
above  lists  will  be  given.  Also  the  background  for  the  research  and 
the  future  perspectives  will  be  discussed. 

BACKGROUND 

The  study  of  soliton  dynamics  in  connection  with  large  Josephson 
tunnel  junctions  has  recently  drawn  considerable  theoretical  [22-27] 
and  experimental  [28,29]  attention*.  Fulton  and  Dynes  [1]  conceived 
the  idea  that  the  Josephson  tunnel  junction  could  support  the  resonant 
propagation  of  a  soliton  (or  fluxon)  trapped  in  the  junction,  the  so¬ 
liton  being  a  2tt  jump  in  the  phase  difference  ( 4> )  across  the  insulat¬ 
ing  barrier  which  separates  the  two  superconductors.  The  moving  soli¬ 
ton  is  accompanied  by  a  voltage  pulse  ( — 4> -^ )  which  can  be  detected  at 
either  end  of  the  junction.  The  dc  manifestation  of  the  motion  is  a 
sequence  of  equidistantly  spaced  branches  in  the  current-voltage  cha¬ 
racteristic  of  the  junction.  These  near-constant  voltage  branches, 
which  were  first  reported  by  Chen,  Finnegan,  and  Langenberg  [30],  are 
known  as  zero-field  steps  (ZFS)  because  they  occur  also  in  the  absence 
of  an  external  maqnetic  field.  In  contrast,  the  so-called  Fiske  steps 
(FS)  are  found  only  when  a  magnetic  field  is  applied  [31], 


* 

References  22-46  are  given  in  the  following  section,  Literature 
cited. 


THE  JOSEPHSON  OSCILLATOR 


An  overlap-geometry  Josephson  tunnel  junction  consists  of  two 
superconductive  metal  layers  (for  example  Nb  and  Pb)  separated  by  a 
thin  insulating  oxide  layer  (NbxOy)  of  uniform  thickness  (t  )  that  is 
small  enough  to  permit  quantum-mechanical  tunnelling  of  electrons.  The 
geometry  is  shown  in  Fig.  1.  The  region  where  the  two  superconducting 
layers  overlap  has  the  length  L  in  the  X-direction  and  the  width  W  in 
the  Y-direction.  Typical  values  are  L  ~  6Aj  and  W  ~  0.8  Aj  where  the 
Josephson  length  A  s  10-1*  m.  Thus  the  overlap  region  is  approximately 
1 -dimensional . 

The  tunnelling  supercurrent  is  described  by  the  two  basic 
Josephson  equations 

j (X,T)  =  jQ  sin<(>  ( 1 ) 

and 


Figure  1.  Josephson  junction  of  overlap  geometry  [32]. 


Here  4>  =  4>{X,T)  is  the  difference  between  the  phases  of  the  order  para¬ 
meter  of  the  two  superconductors/  T  is  laboratory  time,  and  j(X,T) 
is  the  Josephson  current  crossing  the  barrier  per  unit  length  in  the 
X-direction,  j  being  the  maximum  current.  The  voltage  drop  across  the 
insulating  barrier  is  V  =  V(X,T).  Combination  of  (1)  and  (2)  with  Max¬ 
well's  equations  yields  the  following  partial  differential  equation 

<Lp/Rp) 'f’xxT  +  ^XX  "  LP  ^TT  “  GLP 


=  (2TTLpjo/<I>o)  (Sin4>  -  jB/jQ)  (3). 

Here  Lp  and  Rp  are  inductance  and  skin  resistance  per  length  unit  of  the 
oscillator.  (Figure  2  shows  the  equivalent  circuit  diagram  for  the 

oscillator) .  The  capitance  and  the  effective  normal  resistance  per 
length  unit  are  represented  by  C  and  G_1,  respectively.  The  externally 
applied  bias  current  per  length  unit  is  jB,  while 

4>0  =  h/2e  =  2.064  x  10-1  5  Wb  is  the  magnetic  flux  quantum.  Introduc¬ 
tion  of  normalized  coordinates,  x  =  X/Xj  and  t  =  Tw0,  yields  the  per¬ 
turbed  sine-Gordon  equation: 


— -  'i 


Figure  2.  Element  of  lumped  transmission  line  equivalent  circuit 
representing  the  Josephson  oscillator  [32]. 


Here  the  Josephson  length  and  the  Josephson  plasma  frequency  are  given 
by  Xj  =  (4>o/2TTjoLp)  *  and  =  ( 7. ti j Q/ <t,c>C )  ^  r  and  the  coefficients  in  (4) 
are  given  by  a  =  G/w  C,  0  =  u  L  /R  and  y  =  j  /j  .  Typical  values  of 

O  O  r  xr  d  O  ^ 

A  and  u>  are  A_  =  1.6  *  10-4  m  and  w  =  5.8  x  1010s_1  such  that 

J  o  J  o 

the  propagation  velocity  of  electromagnetic  signals  (i.e.  solutions 
of  the  linear  wave  equation,  4>xx  -  <j>  =  0,  corresponding  to  (4))  be¬ 

comes  c  =  AjU0  =  9.3  x  10 6  m/s  in  laboratory  coordinates  for  a  typical 

Josephson  oscillator.  In  the  normalized  coordinates,  x  and  t,  this 
velocity  is  of  course  equal  to  unity. 

At  the  ends  of  the  oscillator,  X  =  0  and  X  =  L,  we  apply  the 
following  boundary  conditions: 

(i)  When  no  external  magnetic  field  is  applied  we  approximate  the 
physical  situation  by  an  open-end  condition,  i.e.  zero  current  on  the 
junction  at  the  ends.  Since  the  current  is  proportional  to  4>v  we  get 

4>x(0,t)  =  4>xU»t)  =  0  (5a) 

in  normalized  units  with  £  =  L/Aj.  Condition  (6)  is  only  an  approxi¬ 
mation  since  coupling  between  the  oscillator  and  the  surrounding  micro- 
wave  circuit  is  neglected.  However,  the  condition  leads  to  good  agree¬ 
ment  between  the  computational  results  obtained  for  this  condition  and 
the  experimental  measurements. 

(ii)  When  an  external  magnetic  field  Hex  in  the  negative  Y-direction 
(see  Fig.  1),  is  present  the  boundary  condition  becomes 

V0,t)  =  4>x(*#t)  =  n  (5b) 

where  n  is  the  normalized  magnetic  field  n  =  (-w/j0*j)Hex  • 

The  initial  conditions  used  for  the  computational  modelling  of 
the  oscillator  are 

<Mx,0)  =  F  (x)  $t  (x,  0)  =  G  (x)  (6) 

where  the  functions  F  and  G  are  chosen  such  that  the  stationary  so- 
liton  dynamic  states  of  the  oscillator  are  obtained  in  the  numerical 
computations  without  too  long  transients.  In  the  following  section 
we  shall  discuss  these  soliton  dynamic  states.  In  principle  the  ini¬ 
tial  conditions  can  be  varied  within  certain  limits  without  any  change 
in  the  resulting  stationary  soliton  dynamic  states.  In  practice  we  often 
use  the  final  values  from  a  previous  numerical  solution  of  a  boundary 
value  problem  (4),  (5)  and  (6)  with  a  slightly  changed  set  of  para¬ 
meters.  The  five  parameters  in  the  model  are  a,  0,  y,  £  and  n. 
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THE  DYNAMICS  OF  THE  JOSEPHSON  OSCILLATOR 


The  classical  sine-Gordon  equation, 

<t>xx“  4>tt  -  sin4>  =  0,  has  2n-kinks  and  anti-kinks 

4>(x,t)  =  4  tan  1lexp(±{x  -  ut  -  xQ)//1  -  u2  ]  (7) 

as  soliton  solutions  [33].  Here  u  is  the  constant  velocity  of  the  sol- 
iton  and  xq  is  the  position  of  the  soliton  at  t  =  0. 

The  perturbed  sine-Gordon  equation  (4)  has  similar  soliton  solu¬ 
tions  in  the  looser  sense.  Each  soliton  carries  a  magnetic  flux  quan¬ 
tum.  The  dynamics  of  these  solitons  is  investigated  by  means  of  per¬ 
turbation  theory  in  Ref.  [34].  As  a  result  a  first-order  differen¬ 
tial  equation  for  the  variable  soliton  velocity  for  a  single  soliton, 
u ( t) ,  is  derived 

^  =  ±  \  ti y  ( 1  "  u2)  3/2  -  au(1  -  u2)  -  -  Bu  .  (8) 

Eq.  (8)  expresses  the  balance  between  energy  input  in  the  system  due 

to  the  Y-term  and  dissipation  due  to  the  loss  terms,  otd>.  and  -84> 

t  Txxt 

The  stationary  velocity,  u^,  is  determined  from  Eq.  (8)  by  letting 
du/dt  =  0  and  solving  the  resulting  equation  with  respect  to  u.  In 
typical  computer  experiments  u(t)  rapidly  adjusts  towards  the  station¬ 
ary  velocity,  u^. 

For  a  finite  junction  with  open-end  boundary  conditions  (5a)  it 
is  easy  to  show  that  solitons  are  reflected  into  antisolitons  and 
vice  versa  at  the  boundaries.  The  bias  current,  y,  drives  the  soliton 
in  the  negative  x-direction  until  it  is  reflected  into  an  antisoliton 
at  x  =  0.  The  antisoliton  is  driven  in  the  positive  x-direction  and 
reflected  into  a  soliton  at  x  =  l  and  a  new  cycle  of  this  stationary 
state  is  initiated.  We  shall  designate  such  a  stationary  state  a  sol¬ 
iton  dynamic  state.  The  periodic  motion  of  the  soliton  on  the  os¬ 
cillator  is  responsible  for  the  emission  of  electromagnetic  radiation, 
typically  in  the  GHz-range,  from  the  oscillator.  Figure  3  shows  a 
computer  picture  of  part  of  the  oscillation  cycle  in  the  soliton  dy¬ 
namic  state  with  one  soliton.  In  the  inset,  <J>  (#,t)  is  shown  as  a  func¬ 
tion  of  t  for  50  time  units.  This  quantity  is  proportional  to  the  vol¬ 
tage  on  the  oscillator  at  the  end  at  x  =  £  according  to  (2) .  The  DC- 
component  of  the  voltage  <|)^(S.,t)  has  been  computed  for  different  va¬ 
lues  of  the  applied  bias  current  y  in  Eq.  (4) .  The  resulting  curve 
shows  agreement  with  experimentally  measured  1st  ZFS  branch  of  the 
IV-characteristic  as  seen  in  Fig.  4.  Similarly  when  two  or  three 
solitons  are  present  on  the  oscillator  the  2nd  ZFS  and  the  3rd  ZFS 
branches  of  the  IV-characteristic  are  obtained. 
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Figure  3.  Computer  solution  of  the  perturbed  sine-Gordon  equation 
(4)  with  a  =  0.05,  0  =  0.02,  y  =  0.35.  Boundary  conditions  (5a)  with 
i  =  6.  Initial  conditions  (6)  with  one  soliton.  The  inset  shows 
4>t(A,t)  [32]  . 

Computational  Fourier  analysis  of  t  » t )  provides  the  power 
spectrum  for  the  radiation  from  the  oscillator.  The  basic  frequency 
is  given  by 

f  =  u/2£  .  (9) 

Also  the  computational  power  spectrum  shows  agreement  with  experimen¬ 
tal  measurements  of  the  power  spectrum  [35-36], 

The  presence  of  the  loss  term-B<)>  in  Eq.  (4)  permits  soliton 
dynamic  states  in  which  two  or  more  s8¥ltons  travel  together  in 
bunches  [34]. 

Figure  3  illustrates  the  2-soliton  case  for  different  parameter 
values.  For  small  values  of  y  (Fig.  3a)  the  two  solitons  travel  in 
a  symmetric  configuration  -  i.e.  soliton  and  antisoliton  in  opposite 
directions.  For  higher  values  of  y  (Fig.  3b)  the  two  solitons  travel 
in  a  bunched  state  -  i.e.  two  solitons  in  the  same  direction  followed 
by  two  antisolitons  in  the  opposite  direction.  There  is  a  gradual 
transition  from  the  symmetric  state  to  the  2-bunch  state  as  the  bias 
current  y  is  increased  and  vice  versa  as  y  is  decreased.  In  Ref.  [37] 
it  has  been  shown  that  the  Hamiltonian  for  two  (undeformed)  solitons 
(on  an  infinite  junction)  has  a  local  minimum  for  a  finite  separation 
between  the  solitons.  This  separation  equals  the  length  of  the  junc¬ 
tion  for  the  value  of  y  at  which  the  transision  between  the  two  sol¬ 
iton  dynamic  states  occurs.  Ref.  [32]  reports  on  the  following  hyste¬ 
resis  phenomenon:  For  increasing  (decreasing)  bias  current  y  the  transi 
tion  from  symmetric  to  bunched  mode  (vice  versa)  occurs  at  smaller 
(higher)  values  of  y. 


Idc(mA) 


Figure  4.  DC-voltage  versus  applied  bias  current  showing  the  first 
three  zero  field  steps.  Circles  indicate  computational  results  and 
solid  lines  represent  experimental  results  [32],  For  the  numerical  re¬ 
sults  we  have  used  a  =  0.05  and  5  =  0.02  in  (4) ,  i  =  6  and  n  =  0  in 
(5),  and  one,  two  and  three  solitons  in  the  initial  conditions  (6) 
on  ZFS  1,  ZFS  2,  and  ZFS  3,  respectively.  Furthermore  y  ~  I.  and 

<<j>t(o,t)>  ~  vdc.  ac 


Fig.  6  shows  the  soliton  trajectories  in  the  xt-plane  correspond¬ 
ing  to  different  soliton  dynamic  states.  The  diagram  ZFS  1  corresponds 
to  Fig.  3.  while  the  diagrams  ZFS  2  (Symmetric)  and  ZFS  2  (2-Bunch) 
correspond  to  Fig.  5a  and  b.  The  diagrams  marked  ZFS  3  (3-Bunch)  and 
ZFS  3 (2-Bunch,  1  Free)  correspond  to  soliton  dynamic  states  found 
experimentally  and  computationally  in  Ref.  [36]  where,  respectively, 
three  solitons  travel  in  one  bunch,  and  two  solitons  travel  together 
and  one  soliton  travels  alone  with  a  slightly  deviating  velocity,  u. 
The  soliton  dynamic  states  corresponding  to  the  other  diagrams  in  Tig. 
will  be  discussed  in  the  following  sections  which  treat  the  results 
obtained  under  the  present  contract  work. 


Figure  5.  Computer  solutions  of  the  perturbed  sine-Gordon  equa¬ 
tion  (3)  with  a  =  0.05  and  8  =  0.2.  Boundary  conditions  (5a)  with 
l  =  6.  Initial  conditions  (6)  with  two  solitons.  (a):  y  =  0.125,  (b) 
y  =  0.3.  The  insets  show  <j»t(Jl,t)  [32], 
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Figure  6.  Soliton  dynamic  states.  ZFS  =  zero  field-step  correspond 
ing  to  boundary  conditions  (5a).  FS  =  Fiske  step  corresponding  to 
boundary  conditions  (5b) . 


LINEWIDTH  OF  ELECTROMAGNETIC  RADIATION  FROM  JOSEPHSON  OSCILLATOR 


In  the  frequency  spectrum  of  the  electromagnetic  radiation  from 
the  Josephson  oscillator  the  width  of  the  lines  is  very  narrow.  An 
experimental  measurement  shows  a  physical  linewidth  as  narrow  as 
1  kHz  of  a  line  at  10  GHz  [38].  The  very  well-defined  frequency  of  the 
radiation  is  a  technologically  important  feature  of  the  Josephson  junc¬ 
tion.  The  narrow  linewidth  is  a  consequence  of  the  coherence  of  the  soli- 
tonic  excitation  of  the  junction. 

In  Ref.  [ 1 ]  of  the  present  contract  soliton  perturbation  theory  is 
used  to  calculate  the  soliton  oscillator  linewidth  arising  from  soliton 
interactions  with  background  radiation.  The  paper  treats  the  line  oscil¬ 
lator  illustrated  in  Fig.  1  and  modelled  by  Eqs.  (4-6)  as  well  as  the 
ring  oscillator  illustrated  in  Fig.  2  of  Ref.ll]'-  In  the  latter  case  the 
boundary  conditions  (4)  are  replaced  by  the  periodic  boundary  conditions 


4>t  ( 0  ,  t)  =  4>tU,t) 
4>x<0,t)  = 


(10) 


where  l  now  denotes  the  circumference  of  the  ring  oscillator.  The  pertur¬ 
bation  analysis  is  based  on  the  ansatz  that  the  soliton  is  given  by 


4>o  =  4 tan-1  exp  [  (x  -  X(t)  )  //l  -  x2  ]  (11) 

where  X  =  X(t)  describes  the  soliton  trajectory  in  the  xt-plane  and 
the  derivative  u  =  X(t)  is  the  velocity.  Following  the  analysis  of 
Ref.  [34]  the  perturbation  theory  is  expanded  to  second  order.  Detailed 
calculations  in  the  case  of  an  oscillator  that  is  long  compared  with  the 
Josephson  length  and  for  which  the  radiation  field  is  thermal  establish 
lower  bounds  for  the  linewidth  of  a  real  oscillator.  These  lower  bounds 
are  not  in  disagreement  with  available  instrument-limited  measurements 
of  X-band  linewidths  less  than  5  kHz.  Fig.  3  and  4  of  Ref.  [1]  show  the 
experimentally  measured  linewidths  as  function  of  the  average  soliton 
velocity  and  the  absolute  temperature  of  the  junction  while  Fig.  5  shows 
the  computational  linewidth  at  temperature  T  =  3  K  and  velocity  u  =  0.8. 

A  Hamiltonian  perturbation  theory  was  developed  for  the  ring  oscil¬ 
lator  in  Ref.  [6],  Here  the  ansatz  for  the  travelling  wave  was  given  by 

4>  =  sin-1  [±cnU,k)  ]  (12) 

where  £  =  (x-ut)/k(1  -  u2)^  and  cn(£,k)  is  a  jacobian  elliptic  function 
with  modulus  k.  4)  given  by  (12)  is  an  exact  solution  to  the  unperturbed 
sine-Gordon  equation  <j>Xx  ”  4>tt  “  sin4>  =  0  with  periodic  boundary  condi¬ 
tions  (10).  The  perturbation  theory  was  developed  for  the  elliptic  func¬ 
tion  in  the  same  manner  as  in  Ref.  [34].  As  a  result  an  ordinary  diffe¬ 
rential  equation  like  (8)  was  obtained.  However,  the  coefficients  y( 
a  and  8  were  replaced  by 
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Y '  =  2yk3/A 

a'  =  2ak2E(k)/A  (13) 

S'  =  26  1(2  -  k2 ) E (k)  -  2(1  -  k2)K(k)]/A 

with 

A  =  2k2E  (k)  +  (1  -  k2 )  JL2  4K  (k)  . 


Here  K(k)  and  E(k)  denote,  the  complete  elliptic  integrals  of  first 
and  second  kind.  The  stationary  velocity  was  determined  as  function  of 
circumferential  length  for  different  values  of  loss  parameters  a  and  6 
and  as  function  of  bias  y.  Comparison  to  direct  numerical  solutions  shows 
that  Hamiltonian  perturbation  theory  can  be  used  to  predict  the  statio¬ 
nary  one-soliton  velocity  of  the  ring  oscillator.  The  results  are  use¬ 
ful  for  the  interpretation  of  experimental  measurements  of  the  I-V  cha¬ 
racteristic  for  this  device.  Furthermore  the  stationary  solutions  are 
useful  as  initial  conditions  for  numerical  simulations  of  the  ring  os¬ 
cillator  under  different  circumstances. 


In  Ref.  [11]  very  detailed  simulation  studies  of  the  dynamics  in 
long  ring  oscillators  under  the  influence  of  external  microwave  radiation 
field  and  internal  thermal  noise  are  presented.  The  former  situation  is 
modelled  by  inclusion  of  a  sinusoidal  driving  term  in  the  perturbed  sine- 
Gordon  equation  (4)  where  ~y  is  replaced  by  y  +  n(x,t)  with  n(x,t)  = 
n ( t)  =  nQsinftt,  ft  being  the  frequency  and  no  the  amplitude.  The  latter 
case  is  modelled  by  letting  n(x,t)  be  Gaussian  white  noise  with  zero 
mean  <n(x,t)>  =  0  and  autocorrelation  function 


Rn(C,T)  =  <n(x,t)n(x  +  C/t  +  T)>  =  o26(c)6(t) 


(14) 


Here  the  variance  of  the  noise  oz  is  connected  with  the  loss  a  and  the 
absolute  temperature  T  through  n 

6 2  =  4^T/4»ojo  X  j  (15) 

where  kg  is  the  Boltzmann  constant.  The  simulation  algorithm  uses  a 
pseudo  spectral  method  making  heavy  use  of  fast  Fourier  transforms  which 
is  well  adapted  to  vector  processors  (CRAY-1-S)  which  gives  a  speed-up 
factor  in  computing  time  of  typically  22  in  comparison  to  conventional 
high-speed  computers,  and  also  provides  results  with  a  relative  accuracy 
of  less  than  10“ 8  which  is  required  to  study  the  very  narrow  linewidth 
of  such  oscillators.  For  the  sinusoidal  driving  term  the  computational 
results  are  compared  to  results  obtained  by  perturbation  theory  based 
on  the  ansatz 

<Mx,t)  =  4>k  (x,  t)  +  4>°°(t) 


(16) 


where  <p  (x,t)  denotes  the  kink  part  and  4>  (t)  denotes  the  background 
part.  In  the  theory  the  background  motion  becomes  an  effective  driving 
term  for  the  momentum  of  the  kink  part.  As  a  result  the  kink  momentum  is 
determinded  as  a  function  of  time  and  thus  the  fluctuations  in  the  kink 
velocity  and  in  the  revolution  frequency.  Since  the  velocity  u  fluc¬ 
tuates  Eq.  (9)  is  replaced  by 


t  .  + 
■  n-1 


n 


udt  =  l 


(17) 


and 


'n-1 


_n 


=  1/T 


n 


where  Tn  is  the  n 1 th  revolution  period.  We  compare  the  standard  devia¬ 
tion  of  the  revolution  frequency 

Of  =  (<  (fn  -  <fn>) 2>J *  for  0.4  <  n  <  2.0. 

The  results  are  shown  in  Fig.  8  of  Ref.  [11].  The  perturbative  kink  model 
predicts  a  resonance  just  below  the  plasma  frequency  ft  =  1 ,  whereas 
the  numerical  simulation  yields  this  peak  at  a  somewhat  lower  frequency. 
The  reason  for  the  discrepancy  is  attributable  to  the  fact  that  we  have 
used  a  linearized  kink  model.  Presumably,  the  use  of  a  higher  order  ex¬ 
pansion  would  yield  a  behavior  analogous  to  that  of  a  soft  nonlinear 
spring,  thus  reducing  the  discrepancy.  It  has  also  been  shown  later  that 
the  difference  in  scale  of  the  standard  deviation  can  be  removed  by  a 
further  refinement  of  the  perturbative  treatment  [39], 

For  the  Gaussian  white  noise  driving  term  the  numerical  simulation 
results  agree  well  with  experimental  results.  The  model  based  on  Hamil¬ 
tonian  perturbation  theory  was  also  able  to  predict  the  qualitative  de¬ 
pendence  of  the  standard  deviation  on  the  bias. 

In  Ref.  [7]  standard  methods  of  stochastic  processes  are  used  to 
study  the  coupling  of  the  sine-Gordon  system  with  a  heat  reservoir.  Both 
phonons  and  solitons  are  found  to  be  thermalized  in  a  way  such  that  the 
phonons  will  have  an  average  energy  at  kB1!  per  mode  while  solitons  will 
have  an  energy  of  J  kgT .  These  results  are  in  agreement  with  those  ob¬ 
tained  by  using  a  statistical-mechanics  approach  for  a  "dilute"  solu¬ 
tion  gas  [ 40 ] . 
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EFFECTS  OF  EXTERNAL  MAGNETIC  FIELDS  ON  SOLITON  DYNAMICS 

In  Ref.  [4]  the  dynamical  benavior  of  solitons  propagating  in  the 
presence  of  an  applied  magnetic  field  on  a  long  overlap-geometry 
Josephson  tunnel  junction  is  investigated.  The  application  of  external 
magnetic  fields  is  essential  for  tuning  the  frequency  of  the  radiation. 
The  magnetic  field  is  modelled  by  boundary  conditions  (5b).  It  is  de¬ 
monstrated  that  the  soliton  dynamic  state  corresponding  to  the  branch 
of  the  IV-characteristic  for  the  oscillator  denoted  the  1st  Fiske  step 
(FS  1)  is  the  following:  The  soliton  travels  in  the  bias-aided  direc¬ 
tion  which  is  the  negative  x-direction  in  the  FS  1  diagram  in  Fig.  6 
and  reacts  with  the  boundary  condition  (5b)  at  x  =  0.  As  a  result  energy 
is  absorbed  from  the  incident  soliton  such  that  the  minimum  energy  for 
the  sine-Gordon  soliton  (=  8  in  normalized  units)  is  no  longer  available. 
Therefore  no  antisoliton  is  reflected.  Instead  reflection  of  plasma  os¬ 
cillations  is  observed.  They  travel  in  the  positive  x-direction  and 
reach  the  boundary  at  x  =  t  where  a  new  soliton  is  created  due  to  the 
energy  input  caused  by  boundary  condition  (5b).  This  constitutes  the 
first  cycle  of  the  stationary  soliton  dynamic  state.  For  the  initial 
conditions  used  in  this  paper  the  2nd  Fiske  step  which  contains  two 
solitons  is  not  found,  while  the  3rd  Fiske  step  (FS  3  in  Fig.  6)  is 
found  to  consist  of  two  solitons  travelling  in  the  negative  x-direc¬ 
tion  and  one  antisoliton  and  plasma  oscillations  travelling  in  the 
positive  x-direction.  Analyses  of  the  corresponding  of  the  three  first 
harmonics  of  the  computational  power  spectra  confirm  this  picture  of 
the  soliton  dynamics  for  the  junction  biased  on  FS  1  and  FS  3. 

Ref.  [3]  demonstrates  that  the  1st  Fiske  step  possesses  a  branched 
structure.  The  major  portion  of  the  step  corresponds  to  a  simply  periodic 
soliton  oscillation  whereas  the  branches  are  characterized  by  subhar¬ 
monic  generation.  Such  period  doublings  are  important  because  they  may 
be  the  first  step  on  the  road  to  chaos.  Indeed  such  chaotic  behavior  is 
found  by  us  for  overlap  Josephson  junctions  when  the  oscillator  is  sub¬ 
jected  to  a  constant  external  magnetic  field  [5].  The  phenomenon  will 
be  discussed  in  the  following  section. 

In  Ref.  (2]  the  soliton  dynamics  for  different  geometries  of  the 
Josephson  tunnel  junction  are  compared.  So  far  we  have  only  been  con¬ 
cerned  with  the  overlap  geometry  illustrated  in  Fig.  1.  In  the  in-line 
geometry  the  bias  current  enters  in  the  direction  parallel  to  the  long 
dimension  (instead  of  perpendicular  to  the  long  dimension)  of  the  junc¬ 
tion  and  is  limited  by  self-screening  effects  to  the  two  ends  of  the 
junction.  In  the  paper  I-V  characteristics  and  microwave  emission  spectra 
are  calculated  for  the  two  geometries  and  shown  to  be  qualitatively  si¬ 
milar,  although  also  quantitative  differences  are  found. 

Ref.  [9]  is  a  detailed  study  of  I-V  structure  and  the  emitted  X-band 
radiation  from  the  overlap-geometry  Josephson  tunnel  junctions  of  inter¬ 
mediate  length  ( i  =  2)  when  external  magnetic  fields  are  applied.  Di¬ 
rect  computational  solutions  of  the  perturbed  sine-Gordon  equation  (4) 
with  inhomogeneous  boundary  conditions  (5b)  are  compared  to  experimental 
measurements.  Furthermore  cavity  mode  analyses,  both  single  mode  and 
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multimode,  are  shown  to  predict  the  simulation  results  quite  accurately. 
This  result  is  important  because  the  application  of  a  superposition  of 
linear  cavity  modes  is  computationally  cheaper  than  direct  numerical  so¬ 
lution  of  the  boundary  value  problem  of  the  partial  differential  equa¬ 
tion.  The  method  works  for  intermediate  length  Josephson  junctions.  For 
longer  junctions  the  soliton  nature  of  the  excitation  becomes  more  im¬ 
portant  and  must  be  included  in  the  expansion  modes. 

Fig.  1  in  Ref.  [9]  gives  some  of  the  main  results  of  the  paper.  For 
small  values  of  the  magnetic  field  n  the  junction  is  biased  on  ZFS  1 
while  larger  values  of  n  leads  to  a  soliton  dynamic  state  corresponding 
to  the  2nd  Fiske  step.  Both  FS  2  symmetric  and  FS  2  2-Bunch  (indicated 
in  Fig.  6)  are  found.  In  the  former  case  solitons  and  plasma  oscilla¬ 
tions  are  travelling  in  opposite  directions  on  the  junction  at  the  same 
time.  In  the  latter  case  two  solitons  are  travelling  in  a  bunch  together 
and  reflected  into  plasma  oscillations.  Thus  we  have  demonstrated  the 
existence  of  the  bunching  phenomenon  [32]  also  in  the  magnetic  case. 

For  n  ~  1.8  Fig.  1  shows  the  bias  value  at  the  top  of  the  1st  zero  field 
step  and  the  bias  value  at  the  bottom  of  this  step  as  function  of  the 
magnetic  field.  The  top  decreases  and  the  bottom  increases  as  the  mag¬ 
netic  field  is  increased.  At  n  s  1.8  ZFS  1  no  longer  exists. 

For  1.8  ^  n  X,  5.5  the  junction  operates  on  FS  2  in  the  symmetric  mode  or 
in  the  bunched  mode  as  illustrated  in  the  insets  of  the  figure.  On  FS  2 
it  is  only  possible  to  detect  the  top  bias  value  of  the  step.  As  indi¬ 
cated  by  the  arrows  in  Fig.  1  the  stationary  values  of  the  bias  current 
depend  on  the  direction  of  the  field  variation.  A  similar  hysteresis  phe 
nomenon  was  found  in  the  non-magnetic  case  [32].  The  figure  also  shows 
the  agreement  between  the  results  obtained  by  direct  numerical  solution 
and  the  results  obtained  by  the  single  mode  theory  due  to  Kulik  [41],  up 
to  n  ~  4.4,  and  the  results  obtained  by  the  multimode  theory  due  to 
Enpuku  et  al.  [42]. 

Ref.  [9]  contains  numerous  computations  of  time  series  for  the  vol¬ 
tage  at  one  end  of  the  junction,  power  spectra,  and  dependence  of  fre¬ 
quency  components  on  bias  and  magnetic  fields  (illustrated  in  Figs.  2-9) 

Experimental  measurements,  shown  in  Fig.  10-12,  agree  well  with  the 
computational  results.  In  particular,  it  was  found  experimentally  that 
the  first  frequency  component  was  missing  in  the  interval  in  the  mag¬ 
netic  field  strength  where  the  junction  operates  on  FS  2.  Here  the  po¬ 
wer  in  the  first  harmonic  was  found  computationally  to  be  very  low. 

SOLITON  AND  CHAOS  EFFECTS  ON  THE  JOSEPHSON  OSCILLATOR 

As  mentioned  earlier  the  little  distorted  coherent  solitonic  exci¬ 
tation  of  the  Josephson  oscillator  is  the  reason  for  the  narrow  line- 
width  of  the  electromagnetic  radiation  from  the  device.  Chaotic  features 
in  the  nonlinear  dynamics  of  the  oscillator  may  change  this  picture.  In 
particular,  a  combination  of  random  thermal  effects  and  chaotic  effects 


may  have  serious  effects  on  the  frequency  spectrum.  Thus  it  is  desirable 
for  the  technological  applications  to  avoid  operation  of  the  oscillator 
in  the  regions  of  parameter  space  (i,  a,  0,  ,  ,  n)  where  these  phenomena 

occur.  In  Ref.  [5]  we  have  found  a  region  where  chaotic  intermittency 
between  soliton  dynamic  states  occurs .  It  is,  however,  only  a  relatively 
narrow  window.  The  junction  considered  has  >-  =  5,  cx  =  0.252,  p  =  0,  and 
the  magnetic  field  is  n  =  0.25.  For  y  =  0.454  the  junction  is  biased  on 
FS  1  while  for  y  =  0.5  the  junction  is  biased  on  FS  2.  For  intermediate 
^-values  (y  =  0.456-0.490)  the  oscillator  was  observed  computationally 
to  switch  between  FS  1  and  FS  2  intermittently,  giving  rise  to  a  spe¬ 
cial  branch  in  the  IV-characteristic  which  we  denoted  FS  1$  as  indicated 
in  Fig.  b.  Long  computer  simulations  of  the  phenomenon  made  a  detailed 
study  of  the  power  spectrum  possible  and  also  demonstrated  that  the 
switching  can  be  treated  probabilistically  as  a  Poisson  process.  Expe¬ 
rimental  measurements  [43]  have  in  fact  perhaps  revealed  FS  1$. 

For  the  study  of  chaos  it  is  necessary  to  possess  analytical  and 
computational  tools  for  the  detection  of  the  phenomenon.  We  have  worked 
with  the  Painlevd  test  [44]  and  developed  software  for  computation  of 
return  maps.  In  Ref.  [12]  we  apply  these  tools  to  a  periodically  driven 
rf  superconducting  quantum  interference  device  (SQUID)  consisting  of  a 
ring  with  a  single  Josephson  junction.  The  system  is  described  by  the 
ordinary  differential  equation 

<P"  +  e<t>'  +  sin;  =  a(,sinuot  -  ;  )  .  (18) 

Here  prime  denotes  differentiation  with  respect  to  time,  c  is  a  loss 
parameter,  wp  is  the  driving  frequency  and  r,  the  amplitude.  The  per¬ 
turbing  term,  n;,  has  the  effect  of  confining  the  chaos.  As  a  result 
an  almost  1 -dimensional  return  map  is  found.  The  very  delicate  struc¬ 
ture  of  the  dynamics  as  well  as  the  existence  of  coexisting  attractors 
are  demonstrated. 

Besides  detection  of  chaos  prediction  of  the  phenomenon  is  impor¬ 
tant.  Available  here  is  the  method  of  Melnikov  integrals  [45].  So  far 
the  method  has  only  been  devel  ped  for  ordinary  differential  equations 
that  possess  a  homoclinic  orbit.  In  Ref.  [14]  we  apply  the  method  to 
the  following  equations  with  linear  and  quadratic  damping  terms,  ta; 1 
and  k  ((I*)  2  , 

<f"  +  sin4>  =  e(r  +  r  sin  w  t  -  a;') 

(19) 

4"  +  k(<)>')  +  sin;  =  p  +  p.|  sin.ipt, 

respectively.  The  differential  operators  on  the  left  hand  sides  possess 
homoclinic  orbits  in  phase  space.  The  right  hand  sides  are  the  pertur¬ 
bative  terms  giving  rise  to  Smale  horseshoe  chaos  when  the  Melnikov  inte¬ 
gral  vanishes.  In  the  case  of  quadratic  damping  the  existence  of  an  exact 
solution  to  ;"  +  k(4>')‘  +  sin;  =  0  makes  it  possible  to  avoid  the  inclu¬ 
sion  of  the  damping  term  into  the  perturbation.  The  Melnikov  prediction 
of  chaos  becomes  correspondingly  more  accurate  as  is  demonstrated  by 
comparison  to  numerical  solutions  of  (19). 
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OTHER  RESEARCH  IN  CONNECTION  WITH  THE  PRESENT  CONTRACT 

In  connection  with  the  physical  understanding  of  the  soliton  pheno¬ 
mena  in  the  sine-Gordon  system  the  mechanical  analogue  consisting  of 
elastically  coupled  pendula  subjected  to  gravity  [46]  plays  an  important 
role.  In  Ref.  [13]  a  mechanical  analogue  for  the  double-sine-Gordon 
equation 

*xx  ”  ♦tt  "  A1sin*  ‘  —  sin  2  =  0  (20) 

is  proposed  and  used  to  analyze  solitary  solutions  for  arbitrary  para¬ 
meter  values.  Eq.  (20)  applies  to  other  condensed  matter  systems. 

For  the  numerical  solution  of  nonlinear  evolution  equations  like 
the  sine-Gordon  equation  it  is  a  fundamental  question  to  perform  the 
discretization  in  an  optimal  manner.  Ref.  [8]  investigates  this  question 
in  the  case  where  the  evolution  equation  is  integrable.  A  geometrical 
approach  is  used  to  obtain  a  discretization  that  preserves  the  inte- 
grability.  As  an  illustrative  example  the  discrete  Burger's  hierarchy  is 
analyzed.  The  possibility  of  extending  this  procedure  to  soliton  equa¬ 
tions  which  are  also  integrable  is  briefly  discussed. 

PERSPECTIVES 


The  work  done  in  the  present  contract  is  being  continued  under  the 
following  main  themes: 

Instabilities  of  the  steps  in  the  I-V  characteristics  for  the 
Josephson  junction.  For  optimal  operation  of  the  Josephson  oscillator 
it  is  important  to  understand  the  instabilities  of  the  soliton  dynamic 
states.  For  example  why  does  the  spatially  uniform  excitation  of  the 
junction  become  unstable  for  values  of  the  bias  current  which  are  lar¬ 
ger  than  a  certain  critical  value?  As  a  result  a  spatial  structure  is 
formed  (one  soliton  on  ZFS  1,  two  solitons  on  ZFS  2,  etc.).  And  why 
does  this  structure  become  unstable  again  at  larger  values  of  the  bias 
current  such  that  the  step  has  a  maximum  height? 

Chaos  and  noise  rise  due  to  thermal  effects.  The  relationship  be¬ 
tween  these  effects  is  of  crucial  importance  for  the  narrow  linewidth 
properties  of  the  Josephson  oscillator. 

Modal  expansions  for  longer  Josephson  junctions.  Computational 
costs  have  been  reduced  considerably  by  using  cavity  mode  expansions 
for  intermediate  length  Josephson  junctions.  Can  a  similar  reduction  of 
simulation  costs  be  achieved  for  long  Josephson  junctions  by  using 
expansions  in  terms  of  nonlinear  soliton  modes  from  the  unperturbed 
system? 


Coupling  problems.  For  the  practical  application  of  Josephson  os¬ 
cillators  in  thin  film  electronic  networks  the  coupling  between  the 
oscillator,  other  oscillators  and  surrounding  microstrips  is  essential. 
The  dynamics  of  the  junction  equipped  with  boundary  conditions  modelling 
such  couplings  is  very  important. 
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Soliton  perturbation  theory  is  used  to  calculate  the  fluxon  oscillator  linewidth  arising 
from  fluxon  interaction  with  background  radiation.  Detailed  calculations  in  the  case  of  an 
oscillator  that  is  long  compared  with  the  Josephson  length  and  for  which  the  radiation  field 
is  thermal  establish  lower  bounds  for  the  linewidth  of  a  real  oscillator.  These  lower  bounds 
are  not  in  disagreement  with  recent,  instrument-limited  measurements  of  df-band  linewidths 
less  than  S  kHz. 


I.  INTRODUCTION 

In  1973  Fulton  and  Dynes  pointed  out  that  the 
“zero-field  steps”  observed  in  the  voltage-current 
characteristics  of  long  Josephson  junctions  could  be 
ascribed  to  oscillatory  behavior  of  internal  fluxons 
(or  magnetic  solitons).1  Subsequent  observations  of 
microwave  radiation2  led  to  the  hope  that  such 
structures  could  play  a  technically  useful  role  as  os¬ 
cillators  into  the  millimeter  wave  range.3  Recently, 
some  long  Josephson  junctions  of  high  quality  were 
fabricated  and  tested  at  the  University  of  Salerno4 
and  sent  to  the  Physikalische  Technische  Ilunde- 
sanstalt  in  Berlin5  and  the  Technical  University  of 
Denmark6-7  for  detailed  measurements  in  the  mi¬ 
crowave  range.  Comparison  of  these  microwave 
measurements  with  numerical  and  analog  computa¬ 
tions  of  fluxon  dynamics  (based  on  a  structurally 
perturbed  version  of  the  sine-Gordon  equation)8 
confirms  that  the  original  idea  of  Fulton  and  Dynes 
is  correct. 

Among  other  results  emerging  from  these  experi¬ 
mental  studies  has  been  the  observation  of  a  supris- 
ingly  narrow  oscillator  linewidth:  less  than  5  kHz 
(the  instrument  limit)  at  a  fundamental  oscillator 
frequency  of  10  GHz.7  Our  aim  in  this  paper  is  to 
present  a  theory  of  fluxon  oscillator  dynamics 
which  allows  us  to  predict  the  linewidth  of  a  long 
Josephson  junction  oscillator. 

Our  approach  is  based  upon  the  description  of  a 
Josephson  transmission  line  as  the  sine-Gordon 
equation  with  structural  perturbations  that 
represent  dissipation  and  input  of  energy.9,10  We 
extend  a  recently  developed  soliton  perturbation 
theory"  to  second  order  in  a  small  parameter  pro¬ 
portional  to  the  structural  perturbations  in  order  to 
calculate  the  effect  of  background  radiation  on  soli¬ 
ton  dynamics.12  This  calculation  allows  us  to  de¬ 
fine  an  “instantaneous  frequency"  which  leads 


directly  to  an  explicit  formula  for  oscillator 
linewidth  as  a  function  of  the  background  radiation 
in  the  junction.  Such  background  radiation  may  be 
generated  in  several  ways:  (i)  electrical  noise  con¬ 
ducted  to  the  oscillator  through  bias  and  output 
leads,  (ii)  radiation  generated  by  spatial  inhomo- 
geneities  of  the  junction,  (iii)  radiation  generated 
during  reflection  of  a  fluxon  from  the  end  of  the 
junction,  and  (iv)  thermal  noise  in  the  cavity  modes 
of  the  junction.  To  obtain  a  lower  bound  on  oscilla¬ 
tor  linewidth,  we  assume  the  radiation  field  to  be 
entirely  thermal  noise.  Under  this  assumption,  and 
with  some  simplifications,  wc  calculate  suitably 
normalized  values  for  linewidth  as  a  function  of 
temperature  and  average  fluxon  velocity.  The  worst 
(i.e.,  largest)  value  of  linewidth  that  we  calculate 
under  these  assumptions  is  less  than  the 
instrument-limited  value  of  5  kHz.7 

Although  the  work  reported  here  is  related  to  re¬ 
cent  studies  of  chaotic  behavior  in  the  sinusoidally 
driven  nonlinear  pendulum  and  sine-Gordon  equa¬ 
tion,13  we  emphasize  that  our  results  depend  upon 
the  assumption  that  the  trajectory  of  the  fluxon  os¬ 
cillation  is  not  trapped  in  a  region  of  phase  space 
that  contains  a  “strange  attractor.”14  The  above- 
mentioned  numerical  studies8  support  this  assump¬ 
tion. 


II.  DESCRIPTION  OF  OSCILLATORS 

Our  analysis  of  fluxon  oscillators  is  based  upon  a 
previously  developed  theoretical  model  for  the 
Josephson  transmission  line9,10  (JTL),  which  is 
briefly  recapitulated  here  for  the  convenience  of  the 
reader.  Figure  1  shows  a  transmission  line 
equivalent  circuit15  for  JTL  in  which  L  is  series  in¬ 
ductance  per  unit  length  (pul)  related  to  supercon¬ 
ducting  surface  currents,  R  is  series  resistance  pul 
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(b) 

FIO.  1.  (a)  Physical  model  of  the  JTL  (not  to  scale), 
(b)  Transmission  line  equivalent  circuit  of  the  JTL. 


related  to  normal  surface  currents,  C  is  shunt  capa¬ 
citance  pul  related  to  electric  Held  in  the  junction,  G 
is  shunt  conductance  pul  related  to  normal  electron 
conduction  across  the  junction,  T  is  an  externally 
imposed  bias  current  pul,  and,  finally, 
•/0sin(2ir4>/4>o)  is  the  Josephson  current  pul  crossing 
the  junction.  Kirchhoffs  equations  for  this  JTL 
model  lead  to  the  following  partial  differential 
equation  for  transverse  voltage  (F): 

&xxt + ♦jtjt  ~  LC&  rr  —  GL  <I>  r 

=70f*sin(2ir<l>/<I>0)  +  rL  ,  (1) 

where  X  and  T  are  laboratory  space  and  time, 
4>0=A/2e  is  the  flux  quantum,  and 

4)2  /  VdT  .  (2) 

Series  inductance  ( L )  and  shunt  capacitance  (C)  are 


related  to  junction  geometry  by 

.  +  d 

L-Ho  w 

(3) 

and 

d 

(4) 

where  Xi  is  "London"  penetration  depth  for  surface 
currents,  W  is  junction  width,  d  is  thickness  of  the 
barrier  region,  c  is  dielectric  permittivity  for  the 


barrier,  and  pQ  (=4irX  10“ 7  H/m)  is  the  magnetic 
susceptibility  of  free  space. 

For  analysis  it  is  convenient  to  normalize  these 


variables  as  follows: 

*  =  2rr4>/4>0,  (5) 

x  =X/Xj  ,  16) 

t  —  T/rj  ,  (7) 

where  kj  is  the  "Josephson"  penetration  length  and 

TjskjVZC  .  (8) 

With  these  normalizations,  velocity  is  measured  in 
units  of 

U0  =  kj/Tj 

=  1/VTC  ,  (9) 

and  (1)  becomes 

.(10) 

where 


cl=GL/tj  ,  (11a) 

02  L/Rtj,  (lib) 

y22irZ.r/4>oA.J  .  (11c) 

With  or,  0,  and  y=0,  (10)  it  recognized  as  the 

sine-Gordon  equation  with  the  exact  soliton  solu¬ 

tion14 

*=4tan-'  exp  ±j~ •  (,2) 


which  represents  the  propagation  of  a  magnetic  flux 
quantum  or  “fluxon"  along  the  junction.  To  make- 
a  fluxon  oscillator,  one  must  design  a  physical  path 
over  which  the  fluxon  can  execute  periodic  motion. 
Two  examples  are  shown  in  Fig.  2.-  In  the  “line  os¬ 
cillator”  [Fig.  2(a)],  a  fluxon  approaches  one  end.  is 
reflected  as  an  anitfluxon  [change  of  sign  in  (12)], 
propagates  to  the  other  end,  and  is  reflected  as  a 
fluxon,  etc.  In  the  “ring  oscillator”  [Fig.  2(b)]  the 
fluxon  proceeds  at  constant  velocity  around  the 
ring.  In  our  calculations,  an  important  parameter  is 
the  total  path,  /,  over  which  the  fluxon  travels  to 
complete  a  cycle  of  oscillation  normalized  to  kj. 
Thus  for  the  line  oscillator  [Fig.  2(a)] 


(13) 


while  for  the  ring  oscillator  [Fig.  2(b)] 
.  2»rR 

h  ' 


(14) 
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lb) 

FIO.  2.  (a)  Line  oscillator.  (b)  Ring  oscillator. 

The  effect  of  the  term  y  in  (10)  is  to  pump  energy 
into  the  fluxon  motion  while  the  a  and  0  terms  dis¬ 
sipate  energy.  In  the  following  section  we  use  soli- 
ton  perturbation  theory  to  calculate  effects  of  these 
terms  on  the  motion. 


III.  OUTLINE  OF  PERTURBATION  APPROACH 

The  approach  to  sine-Gordon  soliton  perturba¬ 
tion  analysis  in  Ref.  11  begins  with  a  nonlinear 
equation 

N?=e/(?) .  (15) 

where  JV?=0  is  a  completely  integrable  (i.e.,  “soli¬ 
ton”)  equation  0=col(0,0,),  where 

col(x,y)s 
and  e  is  a  small  parameter.  Expanding 

one  finds  that 

Nj0mO ,  (17) 

so  is  an  exact  multisoliton  solution  which  de¬ 
pends  upon  certain  constant  parameters  py  (e.g.,  the 


speeds  and  positions  of  the  solitons).  Thus 

•  ‘  (18i 

If  and  arc  secular  (i.e.,  grow  linearly  in  time) 
the  second  and  third  terms  on  the  right-hand  side 
(RHS)  of  (16)  are  useful  only  for  times  of  order  c‘  ‘ 
and  e~l,  respectively.  To  overcome  this  objection 
one  can  allow  order  e  time  variations  in  the  p/  s  of 
4>qso4>i  and  02  satisfy 

L^i=Fi(^0) ,  (19) 

L4i  =  FtW0,$t) ,  (20) 

where  Ft  and  Ft  acquire  extra  terms  because  of  the 
modulations  of  the  p/s,  and  L  is  a  linearization  of 
N  around  0O-  Now  secular  growth  of  0(  and  02  can 
be  avoided  by  requiring  that 

FtLS-JL1),  (21) 

F2l^af) ,  (22) 

where  -/f^(L*)  is  the  discrete  null  space  of  the  ad¬ 
joint  of  L.  From  (21)  and  (22)  one  obtains  ordinary 
differential  equations  (ODE)  for  the  order  e  and  or¬ 
der  c2  variations  in  the  p/s. 

The  strategy  of  our  calculation  is  as  follows.  Or¬ 
der  e  corrections  obtained  from  (21),  are  used  to  cal¬ 
culate  the  effects  of  a,  0,  and  y  terms  in  (10)  on  the 
steady  motion  of  a  JTL  fluxon.  The  radiation  field 
is  then  determined  from  (19).  This  permits  us  to 
evaluate  the  orthogonality  condition  (22)  which 
gives  ODE's  that  determine  the  effects  of  0,  (radia¬ 
tion  field)  on  the  fluxon  motion.  In  our  picture  it  is 
this  interaction  of  the  fluxon  motion  with  the  radia¬ 
tion  field  that  leads  to  an  instantaneous  frequency 
and  therefore  to  a  nonzero  oscillator  linewidth. 

Our  analysis  proceeds  as  follows  (see  Ref.  1 1  for 
details).  The  exact  single  fluxon  solution  (12)  of  the 
unperturbed  sine-Gordon  equation  is  modified  to 

0o=4tan-,[exp(f)]  ,  (23) 

where . 


fsy(f)[x-jr(f)]  .  (24) 

Thus  X(t)  specifies  the  trajectory  of  the  fluxon  and 
y(t)  its  relativistic  contraction.  Two  elements  of 
L*)  are 


(25) 


Ei  s 


(26) 
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The  conditions  FiiFx  and  F|lF2  imply 


(27) 


In  general  we  can  define  a  (time-dependent)  instan¬ 
taneous  frequency  as 

.*  .  ; 

v(t)=v(+u9/l  .  (35) 


Thus,  for  typographical  convenience,  we  define 


The  rms  derivation  of  v(t)  from  its  mean  value  vt  is 


Xau  . 


(28)  Av=|<[v(/)-vc]J>.v),/1.  (36) 


The  time  dependence  of  X  is  divided  into  order  e 
contributions,  calculated  from  (21),  and  order  e2 
contributions,  calculated  from  (22).  Thus 

i-ij+Xj,  (29) 

where 

4 

x  /  "  /(^0(f))scchf dS  .  (30) 

*  —m 

Xt  =  £{\-ui)i'2 

X  j^?sin0o]sechfdf  . 


(31) 


IV.  GENERAL  CALCULATION 
OF  LINEWIDTH 

From  the  results  of  the  previous  section  we  see 
that,  under  steady-state  oscillator  conditions,  the 
fluxon  speed  is 

X  =  ut  +  «,(/)  ,  (32) 

where  ut  is  a  constant  (power  balance)  velocity. 

The  time-varying  component  u„  arises  from  interac¬ 

tion  of  the  fluxon  with  the  radiation  field  and,  from 
<30)  and  (31)  obeys  the  ODE 

u,^Xi-(Xi)„,  (33) 

w  here  <  )„  indicates  a  time  average. 

If  u,=0,  the  fluxon  executes  a  perfectly  periodic 
motion  over  a  path  /  with  frequency 

v(=ut//.  (34) 


We  take  Av  as  a  convenient  measure  of  oscillator 
linewidth.  Since  the  radiation  field  in  (31)  is  not 
periodic  we  take  the  average  in  (36)  over  a  long  time 
as 


Av= 


lim  4;  f  df[v(f)-vt]2 

r-.T-'o 


1/2 


(37) 


Equation  (37),  together  with  (35)  and  (31),  provides 
a  straightforward  procedure  for  calculating  the 
linewidth  of  a  fluxon  oscillator.  We  do  this  for  a 
particular  example  in  the  following  section. 


V.  THERMAL  LINEWIDTH 
OF  A  SINGLE  FLUXON 
OSCILLATOR 


We  now  turn  to  a  practical  question  of  fluxon  os¬ 
cillator  design:  calculation  of  linewidth  when  the 
radiation  field  is  assumed  to  be  in  thermal  equilibri¬ 
um  with  its  environment.  This  calculation  neglects 
other  sources  of  the  radiation  field  (electrical  noise, 
radiation  emitted  from  the  fluxon,  etc.)  thus  it 
should  give  a  lower  bound  for  realizable  linewidths 
and  some  idea  about  how  the  linewidth  depends 
upon  oscillator  parameters  and  temperature.  The 
analysis  is  restricted  to  a  single  fluxon  oscillation  to 
avoid  analytical  difficulties  associated  with  the  phe¬ 
nomena  of  “bunching."* 

We  employ  (31)  where,  from  (10), 


e/(0i)=a 


(38) 


thus  a  is  our  small  parameter  in  the  perturbation 
analysis.  Since  we  are  assuming  that  the  radiation 
field  arises  because  the  linear  modes  of  the  oscilla¬ 
tor  are  in  thermal  equilibrium,  (31)  takes  the  form 


X 2  *  j(  1 — u  *)i/2  f  "[<*/( ) + ( 1  - u  2)( rj<ftl  )Jtanhf  scchf  Jsechf  d£  , 


where  a  is  a  small  parameter  that  measures  the 
structural  perturbation,  and  rj  is  a  small  parameter 
that  measures  the  amplitude  of  the  radiation  field. 

In  the  following  analysis  we  make  two  simplify- 


(39) 

r . . 

ing  assumptions: 

0=0  (40) 


and 
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/» 1.  (41) 

The  first  of  these  is  not  a  serious  restriction  if  one 
assumes  a  somewhat  larger  value  of  a  to  account 
for  dissipation  in  the  P  term  of  (38). 


We  calculate  the  thermal  radiation  field  as  a 
sum  of  individual  photon  inodes  of  a  cavity  which 
contains  a  single  fluxon  moving  with  constant  velo¬ 
city  u.17  Thus 


a 


A*  V' \  —  u* 

y/lv  —  uk„ 


km-uoH 

vTTiP 


cos(kMx  —  coMl )  —  sin(kax  —  oj,r)tanh 


x-ut 
V'l  —  u1 


(42) 


where 


(43a) 


and 


k. 


Inn 
i  ‘ 


(43b) 


Since  the  second  term  in  the  integral  (39)  is  an 
odd  function  of  £  while  the  first  term  is  even,  the 
contribution  of  the  second  term  is  small  except 
when  next.  Then  the  ratio  of  the  first  to  second 
term  is  of  order  /  and,  under  assumption  (41),  we 
can  neglect  the  second  term.  Thus  (39)  takes  the 
form 


V+i.i  +  ^ 


sechf  d£  . 


(44) 


The  term  y/a  in  (44)  merely  contributes  a  constant 
to  Xi  which  is  absorbed  in  the  power  balance  condi¬ 
tion  it  determines  ut.  Thus  it  does  not  enter  into 
our  calculation  of  Av.  _ 

The  component  of  Jfj  that  depends  on  the  radia¬ 
tion  field  is 


I 

3f«=2a^’«cosl(*»“—  w»)*  +  0«]  *  (45) 


where 


u7(l—  kmAa 

4y/lrr  ( kmu-am ) 


Xsech 


vkm 


( 1— uJ),/2 


Xexp[— *rk,(l— M1),/1] 


.  (46) 


To  calculate  the  mode  amplitudes  |<4a),  we  assume 
a  mode  at  frequency  oj.  to  have  the  energy 


Em  = 


ficoM 


exp  J 

k,T  1 

(47) 


Strictly  speaking,  the  relation  between  E„  and  Am 
should  be  calculated  for  a  cavity  containing  a  flux¬ 
on;  however,  from  inequality  (41)  this  relation  is  the 
same  as  that  for  an  empty  cavity.  Thus  in  normal¬ 
ized  units 


[  %rrftu0kRo>K/n<i>\ 

exp 

-1 

k,Trj 

L  +Cuo‘u»+  t 

(48) 


From  (48),  (43),  and  (37)  we  obtain 
Av  «i*(l—w*)i/2 


2v/2<J>o/ 


(Ogk,, 

1 

sech 

Zk,U-u2)l/1 

exp[— *,tril—  u2)l/i] 

2 

n 

exp 

k.Trj 

-1 

k2  „  ,  ,  2eJ0kjl 

L  -f  Cluolu«  -f  ^  (k'U-otJ1 

1/2 


(49) 
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Equation  (49)  is  the  main  result  of  this  paper.  To 
appreciate  the  dependence  of  linewidth  Av  upon  os¬ 
cillator  parameters  and  temperature  that  it  implies, 
we  tum  to  two  examples  of  JTL  that  have  been 
thoroughly  studied.10  Important  parameters  are 
recorded  in  Table  I.  From  these  parameters  we 
have  plotted  in  Figs.  3  and  4  /A v/a  as  a  function  of 
u  for  several  values  of  temperature  and  /.  Since 
these  calculations  are  rather  insensitive  to  l  (see  Fig. 
5)  we  can  assume  Av«  We  see  that  A v(u)  rises 
to  a  maximum  value  at 

u~0.8  (30) 

and,  as  we  expect,  falls  to  zero  in  limits  u  — »0  and 
1.  The  main  difference  between  N25L  and  N53C  is 
in  the  value  for  Xj,  but  this  has  a  relatively  small 
effect  upon  Av.  We  find,  of  course,  that  Av  falls 
with  decreasing  temperature,  but  it  is  interesting  to 
observe  that  the  curves  Av(u)  show  little  change  in 
shape. 

Dueholm  et  al.1  have  reported  an  instrument- 
limited  measurement  that 

A <3  ,  (51) 

where  A  is  the  linewidth  (in  units  of  kHz)  for  a  line 
oscillator  with. 

a=0.01  , 

/=  12 . 

From  our  calculation  the  thermal  linewidth  in  labo¬ 
ratory  units  is  given  by  A v/tj.  Assuming  /  =  12 
[i.e.,  a/Xj  —  6  in  Fig.  2(a)]  we  find  for  N25L  that 
the  maximum  linewidth  is  equal  to  260  Hz  and  for 
N53C,  the  maximum  linewidth  is  equal  to  530  Hz. 
These  results  are  not  inconsistent  with  (51). 

VI.  CONCLUDING  DISCUSSION 

The  main  result  of  this  paper  is  (49)  which  gives 
Av/a  as  a  function  of  the  oscillator  parameters 
where  l  is  the  total  fiuxon  path  length  for  a  cycle  of 
oscillation,  measured  in  units  of  X,  u  is  the  average 


TABLE  I.  Josephson  transmission  lines. 


Parameter 

N25L 

N53C 

Unit 

a 

0.0052 

0.00555 

«* 

2.3X10* 

1.76x10* 

m/t 

L 

2.1x10“* 

2.5x10“* 

H/m 

C 

0.9X10“* 

1.3X10* 

F/m 

*7 

i.27xio-> 

2.63  xl0“4 

m 

A, 

9.7x10“* 

1.9 

A/m 

0.55X10-'° 

1.5X10“" 

% 

FIG.  3.  Normalized  thermal  linewidth  as  a  function 
of  average  fiuxon  velocity  and  absolute  temperature  for 
JTL  No.  N2SL. 


fiuxon  speed  normalized  to  u0  (  ~Xj/tj ),  and  T  is 
the  absolute  temperature,  in  addition  to  the  JTL 
parameters  Xj,  Tj,  C,  L,  and  J0  defined  in  Sec.  II. 

From  (49)  the  rms  deviation  of  the  oscillator 
linewidth  is  equal  to  A  v/tj  Hz,  where  a  measures 
the  shunt  oscillator  losses  (including  loading).  In 
deriving  (49)  the  following  assumptions  have  been 
made: 

(1)  only  a  single  fiuxon  is  present  in  the  cavity, 

(2)  /  » 1, 

(3)  the  background  radiation  field  is  entirely  ther¬ 
mal,  and 

(4)  surface  losses  W$ui  in  (10)]  arc  neglected. 

Thus  our  calculations  give  a  lower  bound  for  the 

linewidth  to  be  found  in  a  real  oscillator.  Addition¬ 
al  contributions  to  oscillator  linewidth  may  arise 


FIG.  4.  Same  as  Fig.  3  for  JTL  no.  N53C. 
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FIG.  5.  Itsv/a  vs  /  for  T= 3  K,  and  u=0. 8. 

from  excess  electrical  noise  and  radiation  from  the 
fluxon  itself. 

Equation  (44)  shows  an  exact  mechanism  for  in¬ 
fluence  of  electrical  noise  on  the  fluxon  motion 
through  stochastic  behavior  of  the  bias  current  y. 
A  line  oscillator  may  have  larger  linewidth  than  a 
corresponding  ring  oscillator  because  the  kink- 
antikink  reflection  that  take  place  in  a  line  oscilla¬ 
tor  generate  an  additional  component  of  radiation11 
that  is  not  present  in  a  ring  oscillator.  Since  we  see 


no  special  difficulties  in  making  ring  oscillators,  we 
suggest  that  they  be  considered  experimentally. 

Finally,  Figs.  3  and  4  show  Av  rising  to  a  max¬ 
imum  value  around  u=0.8  Although  this  result  is 
obtained  for  thermal  linewidth,  we  feel  that  this 
behavior  should  be  found  when  a  more  general  radi¬ 
ation  field  is  present.  An  experimental  check  of 
this  suggestion  should  be  possible  with  instrumental 
resolution  of  linewidth  only  an  order  of  magnitude 
better  than  that  reported  in  Ref.  7. 
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Numerical  integration  of  the  perturbed  sine-Gordon  equation  describing  a  long  overlap-geometry 
Josephson  junction  in  a  magnetic  field  indicates  a  branched  structure  of  the  first  Fiske  step.  The 
major  portion  of  the  step  corresponds  to  a  simply  periodic  fluxon  oscillation  whereas  the  branches 
are  characterized  by  subharmonic  generation. 

PACS  numbers:  74.50.  +  r,  85.25.  +  k,  84.20.Pc,  84.30.Ng 


The  dynamics  of  fluxons  on  long  Josephson  tunnel 
junctions  has  recently  attracted  considerable  interest. 
Fluxon  propagation  has  clearly  been  demonstrated  to  be  re¬ 
sponsible  for  the  appearance  of  zero-field  steps  (ZFS)1-3  and 
indicated  as  being  associated  also  with  the  appearance  of 
Fiske  steps  (FS)4-6  in  the  current-voltage  (/-  V )  characteristics 
of  such  long  junctions.  Moreover,  associated  with  this  prop¬ 
agation  there  is  an  emission  of  microwave  radiation  of  very 
narrow  linewidth  from  the  ends  of  the  junction,7  suggesting 
the  possibility  of  interesting  electronic  applications. 

In  this  letter  we  report  on  detailed  numerical  investiga¬ 
tions  of  a  perturbed  sine-Gordon  model  of  a  long  Josephson 
junction  in  a  magnetic  field.  The  results  that  have  emerged 
give  further  confirmation  to  the  fluxon  propagation  mecha¬ 
nism  as  being  responsible  also  for  the  FS  and,  more  impor¬ 
tantly,  show  a  branching  of  the  first  FS  in  the  /-  V  plane  with 
the  presence  of  subharmonic  generation  on  the  branches. 

The  mathematical  model  studied  is1 

-4>„  -sin 0  =  cr^,  -y.  (la) 

*,(0,r)  =  *„(L.r)  =  i7.  (lb) 

Here,  ^  is  magnetic  flux  normalized  to  h/2e,  x  the  longitudi¬ 
nal  distance  normalized  to  the  Josephson  penetration  depth 
Aj ,  and  t  the  ti  me  normalized  to  the  inverse  of  the  Josephson 
plasma  frequency  <u0.  The  y  term  represents  a  uniformly  dis¬ 
tributed  bias  current  normalized  to  the  maximum  zero-vol- 
,  tage  (Josephson)  current  /0  appropriate  to  an  overlap  geome¬ 
try.*  The  term  in  a  represents  quasiparticle  loss.  The 
constant  17  is  a  normalized  measure  of  the  external  magnetic 
field  which  determines  the  boundary  conditions  [Eq.  (lb)]  at 
the  two  ends  of  the  junction  of  normalized  length  L.  For  this 
study  L  —  5,  a  =  0.252, 17  =  0.75,  and  0  <  y  <  1 .  These  val¬ 
ues  were  chosen  to  be  similar  to  those  used  in  Ref.  5.  Equa¬ 
tions  (la)  and  (lb)  were  integrated  using  the  implicit  finite 
difference  method  described  in  detail  in  Ref.  3,  with  the 
space  and  time  intervals  both  set  to  0.05.  Numerical  accura¬ 
cy  and  stability  were  verified  by  halving  and  doubling  the 
space  and  time  intervals  in  the  computations  and  by  carrying 
out  the  integration  for  long  periods  of  time  (r  ~  1500)  and 
observing  that  all  measurable  characteristics  of  the  solutions 
remained  stationary. 

Two  types  of  initial  conditions  were  used:  (i)  a  “smooth” 


"  Permanent  address;  Istituto  di  Fisks,  UnivcrtiU  di  Salerno,  I-M100  Sa¬ 
lerno,  Italy. 


initial  condition,  defined  as 

^(x,0)  =  F(x,0)  +  G(x),  (2a) 

*,(x.O)  =  F,(x,0),  (2b) 

where 

F(x,t)  =  4  tan-'|exp[(x  -  2.5  +  «/)/(!  -  u,),/J]|,(2c) 

in  which  the  velocity  u  (0  <  u  <  1 )  is  chosen  by  a  power  ba¬ 
lance  calculation  according  to  the  value  of  y,  and  G  (x)  is  a 
static  solution  of  Eqs.  (la)  and  (lb)  in  which 
rj  =  0.75  —  (0,0);  (ii)  a  “tickled”  initial  condition,  similar 

to  (i)  but  with  the  further  superposition  of  a  packet  of  plasma 
oscillations,  defined  as 

//  {x,t )  =  A  cos(kx  —  cut  )exp[  —  (x  —  x0)VW/2],  (3) 

in  which  k  =  1,  co1  =  1  +  k  2,  0.1  <,4 <0.2,  1  <x0<4,  and 
0.1<JF<0.2. 

The  results  presented  refer  to  the  situation  with  the 
junction  biased  on  the  first  FS.  Our  most  significant  result  is 
indicated  in  Fig.  1 ,  which  shows  the  /-  V  form,  on  an  expand¬ 
ed  scale,  of  this  FS.  Here,  and  in  the  following,  voltage  is 
defined  as  ,  which  represents  the  physical  voltage  normal¬ 
ized  to  fnaj'le.  In  particular,  in  addition  to  the  major  por¬ 
tion  of  the  step,  similar  to  that  observed  by  Erne  et  al .,5  we 
observe  two  lateral  branches.  These  branches  are  character¬ 
ized  by  the  generation  of  subharmonics  in  the  radiation  emit- 


FIO.  1 .  Detail  of  current-voltage  form  of  Arst  Filkc  step  for  17 «  0.75.  Rec¬ 
tangles:  major  portion,  without  subharmonic  generation:  triangles:  lateral 
branches,  with  subharmonic  generation. 
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FIG  2  Dynamics  of  lower  subharmonic  branch  of  first  Fiske  step  for 
Y  =  0  54,  r)  =  0.75.  Time  evolution  of  phased,  voltage  4,  (upper  inset),  and 
power  spectrum  (lower  inset),  (a)  Left  junction  end,  (b)  right  junction  end. 


ted  from  the  junction,  whereas  the  major  portion  has  no  such 
subharmonic  generation.  Figures  2-4  indicate  the  detailed 
dynamics  associated  with  the  labeled  points  oh  the  /-  Fchar- 
acteristic  of  Fig.  1.  In  Figs.  2(a)  and  2(b)  we  show  the  time 
evolution  of  the  phase  <f>  and  the  voltage  <f> ,  at  the  two  ends  of 
the  junction  and  the  related  power  spectra  corresponding  to 
a  point  on  the  lower  branch.  The  power  spectra  were  calcu¬ 
lated  as  follows:  Printouts  of  the  time  evolution  of  <ft,  were 
examined  and  a  tentative  maximum  superperiod  (highest- 
order  subharmonic)  established  by  eye.  Time  intervals  of 
At~  500  of  these  waveforms,  containing  exactly  an  integral 
number  of  such  superperiods,  were  then  fast  Fourier  trans¬ 
formed  using  a  simple  rectangular  window.  For  display  pur¬ 
poses,  all  of  the  spectra  have  been  normalized  to  an  arbitrary 
value  of^  f  =  10"  ".  Solutions  similar  to  those  in  Fig.  2  were 
invariably  arrived  at  from  "tickled"  initial  conditions.  In 
this  connection  it  is  worth  mentioning  that  the  existence  of 
the  branches  seems  to  depend  only  upon  the  fact  of  tickling 
and  not  upon  the  precise  mode  in  which  tickling  is  effected. 
This  fact  was  established  by  varying  the  amplitude,  width, 
and  initial  position  of  the  tickling  wave  packet  in  the  type  (ii) 
initial  condition  and  observing  that  the  state  into  which  the 
solution  evolved  did  not  change  with  these  variations.  As  is 
apparent  from  Figs.  2(a)  and  2(b),  the  fundamental  period  of 
oscillation  is  approximately  12  normalized  time  units.  The 


FIG.  3.  Dynamics  of  major  portion  of  first  Fiske  step  for  y  =  0  54, 
t)  =  0.75.  Time  evolution  of  voltage  4,  and  power  spectrum  at  left  junction 
end. 


presence  of  the  asymmetric  fluxon  propagation  mechanism 
proposed  in  Ref.  4  is  confirmed  by  the  phase  jump  of  2ir  per 
fundamental  period,5  and  by  the  asymmetry  between  the 
time  evolution  at  the  left  and  right  ends  (note  that  the  polar¬ 
ity  of  the  magnetic  field  used  here  is  opposite  that  used  in 
Ref.  5,  for  which  right  and  left  are  interchanged  here  with 
respect  to  that  work).  Particularly  apparent  from  the  insets 
of  Fig.  2  is  the  strong  subharmonic  contei.t  of  this  oscilla¬ 
tion.  This  is  shown  in  detail  in  the  power  spectra  of  Figs.  2(a) 
and  2(b),  which  indicate  the  presence  of  third  and  sixth  sub¬ 
harmonics  as  well  as  multiples  of  these.  In  contrast  with  the 
results  of  Fig.  2,  when  the  “smooth”  initial  condition  is  em¬ 
ployed  with  the  same  value  of  the  bias,  the  simply  periodic 
time  evolution  indicated  in  Fig.  3  is  obtained.  Although 
some  subharmonic  content  is  still  present,  it  is  very  much 
reduced  compared  with  that  shown  in  Fig.  2. 

In  a  similar  way,  higher  up  on  the  FS  is  a  second  branch 
of  the  /-  V  curve,  the  power  spectra  of  the  dynamics  of  which 
are  indicated  in  Figs.  4(a)  and  4(b).  In  Fig.  4(a),  once  again, 
subharmonic  generation  is  observed,  this  time,  however, 
with  only  a  dominant  second  subharmonic.  As  before,  the 
subharmonic  branch  of  Fig.  4(a)  evolves  from  the  type  (ii) 
initial  condition,  whereas  the  simply  periodic  solution  of 
Fig  4(b)  evolves  from  the  type  (i)  initial  condition. 

In  physical  terms,  the  major  portion  of  the  FS  corre¬ 
sponds  to  a  situation  in  which  a  fluxon  propagates  in  the 
field-aided  direction  (here,  from  right  to  left),  is  reflected  at 
the  end  as  a  localized  plasmon  moving  in  the  opposite  direc¬ 
tion,  which  in  turn  is  reflected  at  the  other  end  as  a  fluxon 
which  resumes  propagating  exactly  as  before.4 '  Propagation 
on  the  subharmonic  branches  is  almost  the  same,  except  that 
perfect  periodicity  of  the  overall  process  is  resumed  only 
after  a  certain  number  of  complete  back-and-forth  cycles. 
The  reason  for  the  existence  of  such  multiple  solutions  for 
given  parameter  values  is  not  understood,  but  the  fact  is  con¬ 
sistent  with  the  frequently  noted  experimental  observation 
of  fine  structure  on  FS.*'  In  the  experiments,  wave  packets 
may  be  created  by  thermal  fluctuations  or  imperfections  in 
the  junction  thus  giving  rise  to  creation  of  subharmonics. 

Of  considerable  interest  is  the  manner  in  which  the  so¬ 
lutions  evolve  as  the  bias  parameter  y  is  varied  in  small  incre¬ 
ments  during  the  computation.  As  we  move  along  the 
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FIG.  4.  Dynamics  of  first  Fiskc  step.  Power  spectrum  of  at  left  junction 
end  for  y  =  0.55  and  ij  =  0.75.  (a)  Subharmonic  branch,  (b)  major  portion  of 


branches,  the  relative  ratios  of  the  subharmonic  components 
vary;  however,  until  we  reach  the  juncture  point  of  the 
branch  with  the  major  portion  of  the  FS,  motion  along  the 
branches  is  repeatable  and  reversible.  Once  the  juncture 


point  is  reached,  however,  a  simply  periodic  solution  corre¬ 
sponding  to  the  major  portion  of  the  FS  is  maintained,  which 
is  also  repeatable  and  reversible.  The  only  way  to  proceed 
from  the  major  portion  of  the  FS  to  the  branches  is  through 
tickling  the  solution,  which  causes  a  jump  to  the  branch. 
This  may  be  related  to  the  fact,  as  is  evident  from  Fig.  1 ,  that 
the  branches  have  negative  differential  resistance.  When  the 
bias  is  progressively  increased  on  the  lower  branch,  the  junc¬ 
tion  switches  abruptly  from  the  branch  to  the  major  portion 
of  the  Fs  for  ^>0.5407.  The  same  procedure  on  the  upper 
branch  leads  to  a  jump  to  the  third  FS  for  y>  0.5505  On  the 
major  portion  of  the  step,  the  junction  switches  to  the  third 
FS  for  ^>0.550  and  to  a  static  zero-voltage  state  for 
Y  <  0.528.  Finally,  our  results  suggest  the  possible  existence 
of  further  fine  structure  of  the  branches.  This  is  presently 
being  investigated. 
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The  dynamical  behavior  of  fluxons  propagating  in  the  presence  of  an  applied  magnetic 
field  on  an  overlap-geometry  Josephson  tunnel  junction  of  length  5A .j  having  a  McCumber 
flt  — Sir  is  studied  by  numerical  integration  of  the  circuit  equations  of  a  50-section  lumped 
RSJ-type  (resistive  shunted  junction)  model.  Resonant  propagating  configurations  corre¬ 
sponding  to  the  first  and  third  Fiske  steps  are  found.  The  fundamental  frequencies  and 
power  levels  of  the  radiation  emitted  from  one  end  when  the  junction  is  biased  on  the  first 
and  third  Fiske  steps  and  on  the  first  zero-field  step  are  comparable,  but  a  magnetic  field 
renders  the  power  spectra  at  the  two  ends  of  the  junction  different. 


INTRODUCTION 


The  zero-field  steps,  or  dc  current  singularities, 
that  appear  (also)  in  the  absence  of  applied  magnetic 
fields  in  the  dc  current-voltage  characteristics  of 
Josephson  tunnel  junctions  which  are  long  in  one  di¬ 
mension  with  respect  to  the  Josephson  penetration 
length  kj  have  by  now  been  convincingly  shown  to 
be  associated  with  the  resonant  propagation  of  flux¬ 
ons  in  the  junction,1-1  according  to  the  mechanism 
first  proposed  by  Fulton  and  Dynes.4  In  this  picture 
the  first  zero-field  step  (ZFS),  which  has  a  voltage 
asymptote  of  4>of//,  where  4>0  is  the  magnetic  flux 
quantum,  T  is  the  electromagnetic  wave  velocity 
within  the  junction,  and  1  is  the  length  of  the  junc¬ 
tion,  is  due  to  the  propagation  back  and  forth  along 
the  junction  of  a  single  fluxon;  the  second  ZFS, 
whose  voltage  asymptote  is  twice  that  of  the  first,  is 
due  to  two  fluxons,  etc. 

In  addition  to  the  ZFS,  when  a  dc  magnetic  field 
is  applied  in  the  plane  of  a  Josephson  tunnel  junc¬ 
tion  (either  long  or  short),  a  second  set  of  current 
steps,  called  Fiske  steps  (FS),  is  observed  in  the 
current-voltage  characteristic.  Successive  FS  occur 
with  a  voltage-asymptote  spacing  just  half  that  of 
the  ZFS.  For  short  junctions,  the  theory  of  FS 
developed  by  Kulik3  quite  satisfactorily  accounts  for 
experimental  observations;  however,  this  is  no  lon¬ 
ger  the  case  for  long  junctions.6  Extensions  of 
Kulik’s  theory  to  long  junctions,  framed  in  terms  of 
cavity-mode  interactions,  have  been  formulated  by 
various  authors.7-10  These  analyses  are  intended  to 
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be  applicable  to  both  FS  and  ZFS,  and  they  have 
found  varying  degrees  of  success  in  explaining  ex¬ 
perimental  observations. 

The  idea  of  applying  the  fluxon  propagation 
model  to  the  explanation  of  FS  in  long  junctions  was 
first  suggested,  but  rejected  as  not  physically  feasi¬ 
ble,  by  Fulton  and  Dunkleberger.  It  was  repro¬ 
posed,  with  an  argument  for  its  feasibility,  by 
Samuelsen.'1  The  essential  ingredients  of  this  pic¬ 
ture  are  the  observation  that  an  applied  magnetic 
field  renders  the  junction  dynamical  equation  asym¬ 
metric  through  the  boundary  conditions,11  thus 
rendering  wave  propagation  along  the  junction 
asymmetric,  and  the  observation  that  the  average 
junction  voltage  in  a  fluxon  propagation  mode  de¬ 
pends  only  upon  the  time-averaged  number  of  flux¬ 
ons  present,  so  that,  for  example,  the  first  FS  is  con¬ 
sistent  with  a  situation  in  which  a  single  fluxon  is 
present  for  half  the  time.  Later,  a  numerical  simula¬ 
tion  result  which  supported  Samuelsen’s  hypothesis 
was  reported  by  Dueholm  et  al. 14 

The  purpose  of  this  paper  is  to  contribute  to  a 
clarification  of  the  situation  through  a  numerical 
simulation  of  a  long  Josephson  junction  in  a  mag¬ 
netic  field.  The  results  that  emerge  give  further  sup¬ 
port  to  the  Samuelsen  mechanism.  Specific  predic¬ 
tions  of  the  frequencies  and  power  levels  of  the  mi¬ 
crowave  radiation  emitted  by  a  junction  current- 
biased  on  a  FS,  as  compared  with  the  same  quanti¬ 
ties  with  the  junction  biased  on  a  ZFS,  are  consistent 
with  experimental  measurements.  Moreover,  the  re¬ 
sults  suggest  further  experimental  measurements  to 
check  the  proposed  mechanism. 
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JUNCTION  MODEL  AND  COMPUTATION 
TECHNIQUES 


Tj_  dVj_ 
rA  dr 


=  y+A.J(^y_|— 2^y+^y  +  I) 


The  junction  model  and  computational  procedures 
employed  in  this  study  are  essentially  those  used  in  a 
previous  study  of  ZFS.1  An  electrical  representation 
of  the  model  is  shown  in  Fig.  1.  The  spatial  depen¬ 
dence  of  the  problem  is  taken  into  account  by  con¬ 
sidering  a  50-section  lumped-circuit  approximation 
to  the  junction.  A  junction  length  of  5k j  is  assumed 
throughout  this  study.  Dissipation  is  assumed  to  be 
due  only  to  the  linear  resistor  R  of  Fig.  1 .  Accord¬ 
ingly,  we  assume  simply  that  the  powers  radiated 
from  the  left  and  right  ends  of  the  junction  at  a 
given  frequency  are  equal  to  the  sum  of  the  squares 
of  the  Fourier  components  Vlm+  K*jn  at  that  fre¬ 
quency  of  the  end-point  voltages  VL  and  VK ,  respec¬ 
tively. 

An  overlap-geometry  junction  is  assumed,  imply¬ 
ing  that  the  bias  current,  IB  per  section,  may  be  con¬ 
sidered  to  be  uniformly  distributed  along  the  junc¬ 
tion.15  The  effect  of  an  externally  applied  dc  mag¬ 
netic  field  Bn  in  the  plane  of  the  junction  and  per¬ 
pendicular  to  its  long  dimension  is  modeled  by  ap¬ 
plying  the  dc  current  IF  as  shown  in  Fig.  1.  The  re¬ 
lation  between  the  two  quantities  is  Bcl=fi0lF/w, 
where  fiQ  is  the  permeability  of  free  space  and  w  is 
the  width  of  the  junction  in  the  field  direction. 

The  Josephson  element  J  in  Fig.  1  is  characterized 
by  the  adiabatic  Josephson  equations 


7=/0sin<£  , 

d6  2  w 

-?T  =  —v 


(la) 

(lb) 


Combining  Eqs.  (la)  and  (lb)  with  the  circuit  equa¬ 
tions  of  Fig.  1  results  in,  for  section  1, 


^_dV_ , 
ta  dr 


—  y+ M  +  A.}  (^2 — )  —  sin^j  —  V  j 


for  section  j  ( 1  <]  <  50), 


(2a) 


—  sin 4>j-Vj  , 


(2b) 


for  section  50, 


r,  dV, o 

rA  dr 


—  y— M  +  kji$i9  —  <f>so)  —  sin05Q—  VK  , 


(2c) 


and  for  all  sections  ( 1  <j  <,  50), 


Here  V  is  voltage  v  normalized  to  I0R,  r  is  time  nor¬ 
malized  to  r^Qo/AIoR  (the  inverse  of  the  gap-sum 
frequency),  rc=RC,  y  =  IB/I0,  M  =  IF/I0,  and 
Aj  s4>(/2ffL/0.  The  value  of  the  parameter  rf/rA, 
which  is  a  measure  of  the  dissipation,  is  assumed  to 
be  10  throughout  this  study;  the  value  of  the  param¬ 
eter  M,  which  is  a  measure  of  the  applied  magnetic 
field,  is  fixed  at  8. 

Equations  (2a)— (2d)  correspond  to  a  discretized 
version  of  the  perturbed  sine-Gordon  equation,  with 
boundary  conditions  given  by  </>x(0,t)  =  <f>xU,t) a  Bex, 
as  employed  by  many  other  authors  (see,  e.g..  Ref. 
13).  The  time  normalization  employed  in  Eqs. 
(2a)— (2d)  is  the  same  as  in  Ref.  1,  but  it  differs  from 
that  used  by  other  authors,  who  measure  time  in 
units  of  the  inverse  plasma  frequency 
l/(oJ  =  [C<t>o/2irI0)'/2.  It  is  trivial  to  show  that  the 
relation  between  the  two  time  scales  is  given  by 
tUyrA  =  (2rf/irrA)-l/J.  It  is  likewise  trivial  to  show 
that  the  dissipative  parameter  rc  /rA  =  2PC  /n,  where 
Pc  is  the  usual  McCumber  parameter. 

Equations  (2a)— (2d)  were  integrated  by  means  of 
a  fourth-order  Runge-Kutta  routine  using  a  fixed 
time  increment  of  0.01  rA.  The  problem  was  in¬ 
tegrated  as  an  initial-value  problem,  and  two  types 
of  initial  conditions  were  employed: 


FIG.  1.  Lumped  circuit  junction  model. 


FIG.  2.  Time  evolution  of  FS-1  fluxon  propagation  for  y=0.53  and  M  =  t.  Circles:  phase  <j>;  crosses:  normalized  volt¬ 
age  V.  (a)  t=266t4  and  270 r4;  (b)  r=274r4  and  278r4;  (c)  r=282r4  and  286r4  (d)  t=290t4  and  294r4.  Apparent  discon¬ 
tinuities  are  due  to  printer  discretization. 


(1)  When  searching  for  a  first  point  on  a  given  FS, 
the  simple  phase  and  voltage  distributions  shown  in 
Fig.  2(a)  of  Ref.  1  were  used.  These  distributions 
represent  only  a  rough  approximation  to  a  single 
propagating  fluxon,  and  their  use  constitutes  a 
weakness  in  the  calculations,  as  will  be  discussed 
below.  (2)  When  searching  for  the  upper  (lower)  ex¬ 
tremity  (in  bias  current)  of  a  given  FS,  y  was  in¬ 
creased  (decreased)  by  a  small  increment  and  the  ini¬ 
tial  phase  and  voltage  distributions  were  taken  as  the 
final  distributions  corresponding  to  the  previous  bias 
value.  The  integration  was  continued  until  the  oscil¬ 
lation  settled  into  a  steady  state,  defined  operational¬ 
ly  as  in  Ref.  1.  Normally,  this  required  integrating 
to  about  300r4  (about  9—10  complete  oscillation 
periods). 

After  stopping  the  integration,  the  final  two  oscil¬ 
lation  periods  of  the  voltages  at  the  two  ends  of  the 
junction,  VL  and  V*,  were  recorded  on  a  grid  of  ap¬ 
proximately  500  equally  spaced  points.  The  average 
values  of  these  voltages,  ( )  and  O'*),  were 
determined  using  the  fact  that  the  average  voltage 
over  a  period  is  proportional  to  the  phase  difference 


over  that  period,  through  the  Josephson  frequency 
relation,  Eq.  (2d).  Physically,  (V)  must  be  constant 
along  the  entire  junction;  in  practice,  ( VL )  and 
{  VR  )  were  always  equal  to  within  better  than  0.5%. 

The  Fourier  sine  and  cosine  coefficients  of  VL 
and  VK  were  then  obtained  by  trapezoidal  rule  in¬ 
tegration  of  the  basic  definitions  of  these  coeffi¬ 
cients  over  the  grid.  Since  the  integration  was  per¬ 
formed  over  two  oscillation  periods,  the  first  even- 
order  harmonic  was,  in  fact,  the  fundamental,  and 
the  absence  of  odd-order  harmonics  served  as  a  fur¬ 
ther  check  that  the  oscillation  had  indeed  reached  a 
steady  state. 

FLUXON  OSCILLATIONS  AND  FISKE  STEPS 

Figures  2  and  3  show  a  series  of  snapshots  of  the 
phase  and  voltage  distributions  along  the  junction 
over  a  single  period  of  oscillation  for  two  different 
modes  of  propagation.  In  Fig.  2(a),  a  fluxon  is  lo¬ 
cated  near  the  center  of  the  junction  at  r  =  266r4 
and  is  propagating  toward  the  right  (this  is  an  anti- 
fluxon  by  the  definition  used  in  Ref.  1,  but  the  dis- 
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FIG.  3.  Time  evolution  of  FS-3  fluxon  propagation  for  y 
age  V.  (a)  r=272rA  and  276rA;  (b)  r=280rA  and  284rA; 
discontinuities  are  due  to  printer  discretization. 


tinction  is  arbitrary).  Note  that  there  is  a  net  phase 
difference  of  2ir  along  the  junction  and  that  the  volt¬ 
age  is  a  well-defined  pulse  which  is  everywhere  posi¬ 
tive.  At  r=270rA,  the  fluxon  is  approaching  a  re¬ 
flection  at  the  right  end  of  the  junction.  In  Fig.  2(b), 
this  reflection  has  just  taken  place  at  t=274ta,  after 
which  a  clearly  defined  voltage  pulse  begins  propa¬ 
gating  to  the  left.  Note,  however,  that  the  voltage 
waveform  now  goes  both  positive  and  negative,  and 
that  the  net  end-to-end  phase  difference  is  essential¬ 
ly  zero.  These  facts  suggest  that  the  entity  in  ques¬ 
tion  is  a  localized  plasmon,  or  perhaps  a  plasmon- 
breather  oscillation.16  In  Fig.  2(c)  this  entity  contin¬ 
ues  propagating  to  the  left,  and  at  r=  286rA  it  is  ap¬ 
proaching  a  reflection  at  the  left  end.  In  Fig.  2(d) 
this  reflection  has  just  occurred  at  t=290ta, 
whereupon  propagation  resumes  toward  the  right. 
During  the  reflection,  however,  the  net  end-to-end 
phase  difference  has  increased  again  toward  2»r,  so 
that  at  r  =  294rA  the  oscillation  has  completed  al¬ 
most  one  full  cycle.  As  is  evident  from  Fig.  2,  the 
phase  at  any  point  along  the  junction  advances  by  a 
total  of  2ir  during  one  oscillation  period,  as  com¬ 
pared  with  the  value  of  4ir  for  the  oscillation  associ- 
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=0.60  and  M  =  8.  Circles:  phase  d;  crosses:  normalized  volt- 
(c)  r=288rA  and  292ta;  (d)  t=296ta  and  300ta.  Apparent 


ated  with  the  first  ZFS  (cf.  Fig.  2  of  Ref.  1). 

In  Fig.  3  the  situation  is  somewhat  less  “binary" 
inasmuch  as  the  various  propagating  entities  overlap 
one  another,  but  the  overall  picture  is  still  sufficient¬ 
ly  clear.  In  Fig.  3(a),  at  r=272rA  there  is  a  net 
phase  difference  of  2 it  along  the  junction  and  a  volt¬ 
age  peak  near  section  21,  corresponding  to  a  fluxon 
moving  to  the  right,  but  a  second  fluxon  is  rapidly 
entering  from  the  left  end.  At  t=276ta  two  voltage 
pulses  are  clearly  visible,  and  the  net  end-to-end 
phase  difference  is  about  3.5jt,  corresponding  to 
something  less  than  two  fluxons  in  the  junction.  In 
Fig.  3(b)  this  packet  has  moved  to  the  right  at 
t=280ta,  with  the  first  pulse  approaching  a  reflec¬ 
tion  at  the  right  end  of  the  junction  and  the  second 
one  near  section  21.  During  this  portion  of  the  os¬ 
cillation  the  situation  suggests,  approximately,  the 
propagation  of  two  fluxons  in  a  “bunched”  configu¬ 
ration.17, 11  At  t=284ta  the  leading  voltage  pulse 
has  just  emerged  from  the  reflection  and  is  located 
near  section  40,  moving  to  the  left;  whereas  the  trail¬ 
ing  pulse  is  located  near  section  33  and  moving  to 
the  right.  In  Fig.  3(c)  the  two  voltage  pulses  have 
exchanged  positions  at  r— 288ta,  with  the  leading 
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FIG.  4.  Normalized  current-voltage  characteristic  of 
the  junction.  Solid  curve:  magnetic  field  equal  to  zero, 
showing  first  two  zero-field  steps;  circles:  first  and  third 
Fiske  steps  for  M  =  8.  Ycm  is  the  field-reduced  critical 
current  value. 

pulse  now  located  near  section  27,  moving  to  the 
left,  and  the  trailing  pulse  just  undergoing  a  reflec¬ 
tion  from  the  right  end.  Note  that  the  net  phase 
difference  along  the  junction  is  at  this  point  essen¬ 
tially  zero.  At  r=292r4  the  two  voltage  pulses  con¬ 
tinue  propagating  to  the  left,  and  the  net  end-to-end 
phase  difference  has  become  approximately  ir.  In¬ 
terpretation  of  this  situation  is  not  completely 
unambiguous,  but  the  propagation  of  a  fluxon- 
plasmon,  fluxon-breather,  or  fluxon-plasmon- 
breather  combination  is  suggested.  In  Fig.  3(d)  the 
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FIG.  3.  End-point  voltage  waveforms  of  FS-1  fluxon 
oscillation  for  y**0.S3  and  M  «8. 


Tlht  (t4l 

FIG.  6.  End-point  voltage  waveforms  of  FS-3  fluxon 
oscillation  for  y=0.60  and  M  =  8. 


leading  voltage  pulse  is  just  undergoing  a  reflection 
from  the  left  end  of  the  junction  at  t=296t4,  and 
the  trailing  pulse  is  located  near  section  19  and  still 
moving  to  the  left.  At  r=300r4  the  leading  pulse, 
barely  distinguishable,  is  located  near  section  17, 
moving  again  toward  the  right,  whereas  the  trailing 
pulse  is  approaching  a  reflection  at  the  left  end.  The 
net  end-to-end  phase  difference  is  now  approximate¬ 
ly  l.Srr,  and  the  oscillation  has  completed  almost 
one  full  cycle.  Figure  3  shows  clearly  that  during 
one  full  period  of  this  oscillation  the  phase  at  any 
point  along  the  junction  advances  by  a  total  of  6ir. 

Figures  2  and  3  show  each  of  the  two  modes  of 
propagation  for  a  single  value  of  the  bias  current. 
Each  of  these  two  modes,  however,  exists  over  a  cer¬ 
tain  range  of  y.  This  fact  is  displayed  in  Fig.  4, 
which  shows  the  two  loci  of  points,  indicated  as  cir¬ 
cles,  labeled,  respectively,  as  FS-1  and  FS-3  in  the 
y—{V)  plane.  The  solid  curve,  for  reference,  is  the 
current-voltage  characteristic  of  the  junction  in  the 
absence  of  magnetic  field,  taken  from  Ref.  1.  The 
positions  of  the  two  groups  of  circles  relative  to  the 
first  two  zero-field  steps,  ZFS-1  and  ZFS-2,  confirm 
that  FS-l  and  FS-3  are,  respectively,' the  first  and 
third  Fiske  steps.  The  point  yCM  *t  ( F)  =0,  in¬ 
cidentally,  is  the  value  to  which  the  critical  current 
has  been  reduced  by  the  magnetic  field  M  —  8. 

At  this  point  we  note  that  there  should  also  exist 
a  second  FS  corresponding  approximately  in  voltage 
position  to  ZFS-1,  but  we  have  not  been  able  to  find 
this  FS  numerically.  We  believe  that  the  reason  for 
this  is  that  the  initial  condition  employed  is  not  real¬ 
ly  appropriate  to  finding  FS-2.  We  have  observed  in 
calculating  FS-1  and  FS-3  that  the  outcome  of  the 
simulation  depends,  at  least  sometimes,  quite  strong¬ 
ly  on  the  initial  condition  employed:  If  the  initial 
condition  is  “too  far”  (by  some  measure)  from  the 
final  propagating  configuration,  the  junction  tends 
to  switch  to  some  other  mode  of  propagation. 
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NORMALIZED  BIAS  CURRENT  Y 

FIG.  7.  Power  spectra  (First  three  harmonics)  of  the  ra¬ 
diation  emitted  from  the  junction  ends,  normalized  to 
VL,VK=0.04/ir.  The  frequency  indicated  is  that  of  the 
fundamental  component.  The  junction  is  biased  on  'he 
First  Fiske  step. 


POWER  SPECTRA 

Figures  2  and  3  clearly  show  that  the  two  modes 
of  propagation  associated  with  FS-1  and  FS-3  are 
asymmetric.  This  fact  is  underlined  in  Figs.  5  and 
6,  which  show  the  temporal  evolution  of  the  two 
end-point  voltages  VL  and  VK  in  the  two  cases  cor¬ 
responding  to  Figs.  2  and  3,  respectively. 

Since  we  assume  that  the  power  radiated  from  an 
end  of  the  junction  is  proportional  to  the  square  of 
the  voltage  at  that  end,  it  is  clear  that  Figs.  5  and  6 
imply  that  the  spectra  of  the  radiation  emitted  from 
the  two  ends  are  quite  different.  That  this  is  indeed 
the  case  is  indicated  in  Figs.  7  and  8,  which  show 
the  power  spectra,  calculated  as  described  above,  for 
the  first  three  Fourier  components  of  the  voltages 
VL  and  VK,  corresponding,  respectively,  to  FS-1  and 
FS-3  as  a  function  of  the  bias  current.  The  frequen¬ 
cy  shown  in  the  lower  part  of  the  figures  is  that  of 
the  fundamental  (labeled  1)  Fourier  component  in 
each  case.  The  power  levels  have  been  normalized 
to  an  arbitrary  value  of  VL,VR  =0.04/ir. 

The  main  conclusions  to  be  drawn  from  Figs.  7 


NORMALIZED  BIAS  CURRENT  Y 

FIG.  8.  Power  spectra  (First  three  harmonics)  of  the  ra¬ 
diation  emitted  from  the  junction  ends,  normalized  to 
Fi.F*  =0.04/ir.  The  frequency  indicated  is  that  of  the 
fundamental  component.  The  junction  is  biased  on  the 
third  Fiske  step. 


and  8  are  the  following:  (1)  The  fundamental  fre¬ 
quencies  associated  with  FS-1  and  FS-3  are  the 
same;  moreover,  they  are  the  same  as  that  associated 
with  ZFS-1  [cf.  Fig.  5(a)  of  Ref.  1].  This  fact  has 
already  been  confirmed  by  experimental  observa¬ 
tions.14  (2)  The  power  levels  of  the  radiation  emitted 
by  a  junction  biased  on  a  FS  are  comparable,  at  least 
in  the  field-aided  direction,  with  those  obtained  with 
the  junction  biased  on  a  ZFS  (cf.  Fig.  5  of  Ref.  1 ). 
This  fact  is  also  consistent  with  experimental  obser¬ 
vations.’4  (3)  The  asymmetries  in  the  power  spectra 
of  the  radiation  emitted  from  the  two  ends  of  the 
junction  may  well  be  large  enough  to  be  measurable. 
One  way  of  doing  this  would  be  to  couple  to  a  single 
end  of  the  junction  and  measure  the  radiation  with 
both  polarities  of  the  magnetic  field. 


COMMENTS 

The  present  work  is  intended  as  a  contribution  to 
the  understanding  of  the  dynamics  of  long  Joseph- 
son  tunnel  junctions.  We  believe  that  fluxon  propa- 
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gation  is  the  basic  physical  mechanism  underlying 
all  of  the  various  phenomena  observed  in  connection 
with  such  long  junctions,  and  that  cavity-mode  in¬ 
teraction  analyses  are  essentially  mathematical  ap¬ 
proximations  that,  as  they  become  more  refined, 
tend  toward  the  same  picture,  being  particularly 
suitable  in  the  limit  of  short  junctions.  We  propose 
that  a  detailed  comparison  of  the  results  of  experi¬ 
mental  measurements  on  real  junctions  with  those 
from  careful  numerical  simulations,  similar  to  those 
recently  reported  by  Lomdahl  el  al. 3  in  connection 
with  ZFS,  should  be  an  important  objective  for  fu¬ 
ture  studies. 
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Chaotic  intermittency  between  soliton  dynamic  states  has  been  found  in  a  perturbed 
sine-Gordon  system  in  the  absence  of  an  external  ac  driving  term.  The  system  is  a  mod¬ 
el  of  a  tong  Josephson  oscillator  with  constant  loss  and  bias  current  in  an  external  mag¬ 
netic  field.  The  results  predict  the  existence  of  a  current  step  between  the  first  two 
Fiske  steps  in  the  current -voltage  characteristic.  A  simple  probability  model  demon¬ 
strates  an  asymmetry  in  the  statistical  nature  of  the  switching  in  the  two  directions. 

PACS  numbers:  74.50.  +  r,  05.40. +  J,  84.30.Ng,  85.25. +  k 


Chaos  phenomena  have  been  found  for  the  rf- 
and  dc-current  driven  small  Josephson  junction 
described  by  the  resistively  shunted  junction 
(RSJ)  model.1  Recently  Ben- Jacob  et  al.  and  Yeh 
and  Kao  reported  on  intermittent  chaos  in  the 
numerical  solution  of  this  model.2  For  a  long 
Josephson  junction,  the  perturbed  sine-Gordon 
equation  (SGE)  with  spatially  uniform  or  nonuni¬ 
form  oscillating  driving  forces  and  linear  damp¬ 
ing  also  gives  rise  to  chaos  in  time  and  space- 
time.3  .  Detailed  numerical  investigations  have 
revealed  subharmonic  generation  caused  by  soli¬ 
ton  motion  in  a  long  Josephson  junction  in  a  con¬ 
stant  external  magnetic  field  modeled  by  the 
perturbed  SGE  without  an  external  ac  driving 
term  in  the  current  bias.4  In  the  present  Letter 
we  demonstrate  a  new  chaotic  intermittency 
phenomenon  between  two  dynamic  states  of  this 
model.  The  two  states  correspond  physically  to 
the  first  two  Fiske  steps  (FS1  and  FS2,  respec¬ 
tively)  in  the  current-voltage  characteristic  of 
the  junction. 

The  mathematical  model  studied  is5 

-  sin<p  =  a<p, -y,  (la) 

vM(o,t)-  ?,(/.,  0  =  n.  (lb) 

Here  ip  is  the  usual  Josephson  phase  variable, 
x  is  distance  normalized  to  the  Josephson  pene¬ 
tration  depth  A  j,  and  t  is  time  normalized  to  the 
inverse  of  the  Josephson  plasma  frequency 
The  y  term  represents  a  uniformly  distributed 
constant  bias  current  normalized  to  the  maximum 
zero-voltage  (Josephson)  current.  The  term  in 
a  represents  quasiparticle  loss.  The  constant  ij 
is  a  normalized  measure  of  the  external  magnet¬ 
ic  field  which  determines  the  boundary  conditions 
(lb)  at  the  ends  of  the  junction  of  normalized 
length  L.  In  this  study  L  =  5,  a  =  0.252,  rj  =  1.25, 
and  y  is  varied  in  the  range  y  =  0.45  -  0.55.  Equa¬ 
tions  (1)  were  integrated  from  appropriate  Initial 
conditions  Isimilar  to  Eqs.  (2)  of  Ref.  4  or  con¬ 


tinuations  from  runs  done  at  nearby  points  in 
parameter  space]  for  long  periods  of  time  (typi¬ 
cally  t  ~  10000)  by  means  of  the  implicit  finite- 
difference  method  described  in  detail  in  Ref.  5, 
with  space  and  time  intervals  set  equal  to  0.05 
and  0.025,  respectively.  Numerical  accuracy 
and  stability  were  verified  by  halving  the  space 
and  time  intervals. 

In  a  narrow  range  of  relatively  low  y  values  (y 
=  0.450-0.454)  the  solution  develops  stably  into 
the  FS1  solution  illustrated  in  Fig.  1(a).  In  physi¬ 
cal  terms  FS1  corresponds  to  a  situation  in  which 
a  soliton  propagates  in  the  field-aided  direction, 
and  is  reflected  at  x  =0  as  a  localized  plasma 
wave  because  of  the  energy  loss  at  this  boundary. 
This  plasma  wave  then  moves  in  the  opposite  di¬ 
rection  and  is  reflected  as  a  soliton  at  x  =  L  as  a 
result  of  the  energy  injection  here  by  the  magnet¬ 
ic  field.  [  An  unambiguous  interpretation  of  Fig. 
1(a)  in  the  schematic  terms  shown  in  Fig.  1(c) 
may  be  established  by  examining  together  <y(x,  t ), 
yx(x,t),  and  y,(x,  t).  |  The  soliton  then  resumes 
propagation  as  before  the  ref  lections. ■‘,6, 7  On  the 
average  this  cycle  lasts  the  time  Tfsi-12  (at  y 
=  0.454).  For  y  <0.450  the  junction  switches  into 
a  static  zero-voltage  state  where  <p,  =  0. 

For  relatively  high  y  values  (y  =  0.500  -  0.540) 
the  solution  develops  stably  into  the  FS2  solution 
corresponding  to  a  situation  in  which  a  soliton 
and  a  localized  plasma  wave  travel  in  opposite 
directions  at  the  same  time  [  shown  numerically 
in  Fig.  1(b)  and  schematically  in  Fig.  1(c)].  The 
average  length  of  this  cycle  was  estimated  to  be 
ffS2 a  6  (at  y  =  0.500).  For  y  >  0.540  the  junction 
switches  to  FS3. 

For  intermediate  y  values  (y  =  0.456-0.490) 
the  junction  exhibits  chaotic  intermittency  be¬ 
tween  FS1  and  FS2.  The  intermittency  is  shown 
in  Fig.  2  for  y  =  0.480,  where  it  is  evidenced  by 
changes  in  the  average  slope  of  p  versus  f  (on 
the  FS1  portions  (*?,)- 0. 6  and  on  the  FS2  por- 
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FIG.  1.  <px(x,t)  on  (a)  FSl  and  (b)  FS2  portions  of 
solution  In  Intermlttency  region  (y  =0.480).  (c)  Sche¬ 
matic  representations  of  (a)  and  (b)  showing  sollton 
(solid  lines)  and  plasma-wave  (dashed  lines)  trajec¬ 
tories.  A  switch  from  FS1  to  FS2  Is  seen  at  t  a  1625. 


FIG.  2.  Phase  vs  time  at  the  left  end  ( x  =0)  showing 
Intermlttency  between  FSl  and  FS2  at  y  =0.480.  Inset: 
Power  spectrum  of  (0,  /),  normalized  to  tp,1  ■  4 
x  10*T,  for  t  =1600-3238,  with  frequency  measured 
In  units  of  the  Josephson  plasma  frequency. 

on  FSl  and  FS2,  respectively.  The  average  times 
the  junction  remains  on  FSl  and  FS2  were  found 
to  be  rFSl  —46  and  fFS,  352  (at  y  =  0.480).  We 
have  found  that  7*FS,  and  rFsa  remain  stationary  in 
time  at  these  values.  The  corresponding  lines 
are  barely  visible  in  the  low  end  of  the  power 
spectrum. 

For  a  particularly  long  time  interval  shown  in 
Fig.  3  (f  =  1902-2079)  the  junction  remained  on 
FS2.  Around  f  =  1960  the  development  of  aperio- 
dicity  is  observed.  The  insets  (a)  and  (b)  show 
the  power  spectrum  for  the  time  intervals  1902- 
1951  and  1984-2079,  respectively,  again  ob- 


/ 


tions  (<P,)*1.0).  The  switch  from  FS2  to  FSl  oc¬ 
curs  when  the  plasma  oscillation  fails  to  generate 
a  sollton  at  x  -  L.  Conversely,  a  switch  from 
FSl  to  FS2  may  occur  when  an  extra  soliton  is 
generated  at  x  =  L  during  the  FSl  cycle.  The 
power  spectrum  (inset  in  Fig.  2)  of  <Pt( 0,  f),  nor¬ 
malized  to  <P,a  =  4x  10' 7 ,  was  obtained  by  a  fast 
Fourier  transform  over  the  time  interval  1600- 
3238  by  use  of  a  Hamming  window.8  Since  loading 
effects  are  not  included  in  the  model,  the  power 
levels  in  this  spectrum  represent  ideal,  available 
values.  Translation  into  physically  measurable 
quantities  requires  a  knowledge  of  the  junction- 
to-mlcrowave  circuit  coupling  (to  set  the  physi¬ 
cal  power  scale)  and  the  junction  plasma  fre¬ 
quency  (to  set  the  physical  frequency  scale).  The 
spectrum  shows  numerous  frequency  components 
among  which  the  two  dominant  lines  at  /*'  =  0,10 
and  0.17  may  be  ascribed  to  the  soliton  motion 
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FIG.  3.  Phase  vs  time  at  the  left  end  (x  =0)  for  a 
particularly  long  operation  on  FS2  at  y  =0.480.  Insets: 
Power  spectra  of  normalized  to  10*T, 

for  (a)  /  =  1902-1951  and  (b)  /  =  1984 -2070.  The  fun¬ 
damental  frequency  Is  /#  =  0.17. 
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FIG.  4.  Abscissa:  Length  of  time  Interval  on  FS1 
and  FS2  measured  In  terms  of  number  of  cycles,  n,  of 
lengths  rFS,  and  rFSi,  respectively.  Ordinate:  Dots 
show  number  of  Intervals,  N((*),  with  <  =  1,2,  shorter 
than  or  equal  to  it.  Full  curves  result  from  theoretical 
estimates  with  pn(t)  =  (1.0  x  10'5)f,  pnU)  =  (1.2  x  10‘4)< 
for  it  <  6,  p2l(t)  =  (0.95  x  10”*)<  for  it  >  6.  Bias  y  =0.480. 
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FIG.  5.  Detail  of  current-voltage  characteristic  for 
L  =  5,  or  =  0.252,  t|  =  1 .25.  FS1  (circles),  FSlJ  (squares), 
and  FS2  (triangles).  Error  limits  for  (</>,)  are  indicated 
by  horizontal  bars  when  larger  than  0.01. 


tained  by  use  of  a  Hamming  window.  In  inset  (a) 
we  observe  the  dominant  frequency  at  f*‘s  0.17 
and  a  relatively  low  content  of  subharmonics.  In 
inset  (b)  subharmonics  have  developed  at  approx¬ 
imately  3  and  f  of  this  frequency.  A  similar 
building  up  of  subharmonics  is  not  seen  while  the 
junction  operates  on  FS1. 

Figure  4  (dots)  shows  the  accumulated  distribu¬ 
tions,  Afj(w)  and  N2M,  of  the  lengths  of  the  time 
intervals  the  junction  operates  on  FS1  and  FS2 
for  217  switches  between  the  two  steps  during  a 
run  of  over  10000  time  units.  Here  n  is  f/fFS, 
and  f/rFS2,  respectively.  To  analyze  this  situa¬ 
tion,  let  the  probability  that  the  junction  switches 
from  FS1  to  FS2  in  the  time  interval  [t  ,t  +d/], 
where  t  is  the  time  after  the  last  switch,  be 
puit)dl.  Furthermore,  if  all  the  switches  are  in¬ 
dependent,  then  the  switching  is  a  Poisson  proc¬ 
ess.  Consequently,  the  probability  Pytf)  that  the 
oscillator  switches  from  FS1  to  FS2  within  the 
time  t  becomes  P,(f)=  1 -exp[-/0'/>12(T)dT].  As¬ 
suming  pu(J)  =pt ,  with p  =  l.Ox  10'3,  we  get  for 
FS1  the  fit  Nt-Pi(nlfSl)Nlm„,  where  Nlmix  is  the 
total  number  of  intervals  on  FS1,  shown  in  Fig. 

4  (full  curve).  The  agreement  between  the  nu¬ 
merical  data  and  this  simple,  but  arbitrary, 
probability  model  is  quite  good. 

For  FS2  a  jump  of  63  in  N2(n)  at  n  =  6  is  ob¬ 
served  in  our  10000-time-unit  run.  Only  by  use 
of  a  transition  probability,  p2l(t),  containing  0.96 
times  a  unit  impulse  at  «  =  6,  besides  a  linear 
term,  do  we  obtain  the  fit  shown  in  Fig.  4.  The 
difference  between  the  two  results  clearly  dem¬ 
onstrates  an  asymmetry  in  the  statistical  nature 
of  the  switching  in  the  two  directions.  The  rea¬ 
son  for  this  fact  is  not  known.  A  possible  cause 


may  be  an  interference  between  the  subharmon¬ 
ics  y  and  §  which  build  up  on  FS2,  as  seen  in  Fig. 
3,  but  the  question  certainly  requires  further 
study. 

Finally,  Fig.  5  shows  the  resulting  current- 
voltage  characteristic  for  the  y  interval  covering 
FS1,  the  intermittency  region,  and  FS2.  At  the 
average  voltage  (  )  a0.83  a  jump  in  the  current 

(from  y  =  0.462  to  y  =  0.480)  occurs.  We  pro¬ 
pose  the  name  “FSl£"  for  this  portion  of  the  char 
acteristic.  We  have  checked  that  the  values  of 
(<Pt)  are  stationary  in  time.  Thus  it  should  be 
possible  to  detect  FSl|  experimentally.  Recent 
measurements  by  Cirillo,  Costabile,  and  Par- 
mentier8  have  in  fact  perhaps  revealed  such  struc 
tures. 
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A  hamiltonian  perturbation  theory  is  developed  for  the  perturbed  sinc-Gordon  equation  with  periodic  boundary  con¬ 
ditions  modelling  the  Josephson  ring  oscillator.  Stationary  fluxon  velocities  arc  determined  as  function  of  length,  loss  and 
bias  parameters. 


The  circular  Josephson  oscillator  was  originally 
proposed  by  Scott  and  McLaughlin  [1  ]  who  pointed 
out  that  this  oscillation  structure  may  be  of  technical 
importance  at  electromagnetic  wavelengths  of  100 
pm  or  less.  This  oscillator  structure  has  in  fact  been 
included  in  the  current  programmes  of  several  ex¬ 
perimental  groups  [2] .  The  mathematical  model  for 
the  circular  Josephson  oscillator  with  losses  and  bias 
current  is  a  perturbed  sine-Gordon  equation  J3J  with 
periodic  boundary  conditions  and  fixed  winding  num¬ 
ber.  These  boundary  conditions  close  the  circular  os¬ 
cillator  in  a  perfect  matching  and  thus  permit  undis¬ 
turbed  fluxon  motion  on  the  oscillator.  Furthermore 
the  boundary  conditions  are  ideal  for  spectral  method 
numerical  studies  of  the  radiation  from  the  oscilla¬ 
tor.  However,  perturbation  methods  for  fluxon  dy¬ 
namics  so  far  have  only  been  used  for  infinitely  long 
Josephson  junctions  [1  ].  In  the  present  note  41  the 
perturbation  method  is  extended  to  the  finite  case 
with  periodic  boundary  conditions. 

The  normalized  perturbed  sine-Gordon  equation 
[3]  with  periodic  boundary  conditions  can  be  written 

*  Supported  by  the  Danish  Council  for  Scientific  and  In¬ 
dustrial  Research  and  by  United  States  Army  through  its 
European  Research  Office. 

41  Based  on  a  master's  thesis  by  one  of  the  authors  (I 


<t>x(0,t)  =  <t>x(U),  0,(0,  o  «0, (/./),  (1) 

where  the  a  term  represents  quasi-particle  loss  across 
the  barrier  and  the  0  term  the  surface  impedance  of 
the  superconductor.  The  7  term  is  the  bias  current. 
The  circumference  of  the  circular  transmission  line 
normalized  to  the  Josephson  length  is  denoted  /  in 
the  periodic  boundary  conditions.  According  to 
hamiltonian  perturbation  theory  |1J  the  hamiltonian 
for  the  unperturbed  sineGordon  equation, 

/ 

11  =  f  +  j0,2  +  1  -  COS0)djr, 

0 

satisfies  the  differential  equation 
dll  ' 

-57-  =  ~  /  (ft0,2  +  +  70,)dx  (2) 

0 

for  small  values  of  a,  0,  and  7. 

The  travelling  wave  solution  to  the  unperturbed 
sineGordon  equation  is  given  by  [4] 

0  =  2  sin-1  [±cn(£,fc)] ,  (3) 

with  {  =  (x  -  ut)/k(\  -  u2)'/2.  Here  cn(£,A;)  is  a 
jacobian  elliptic  function  (5)  with  modulus  *.  Plus 
and  minus  sign  refer,  respectively,  to  fluxons  and 
antifluxons.  The  velocity  of  the  wave  is  denoted  u. 

If  the  modulus  satisfies  the  condition 
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11(1  -u2)'12  =  2kK(k),  (4) 

the  periodic  function  <j>  gets  the  period  /,  and  <t>  rep¬ 
resents  a  single  fluxon  (or  antifiuxon).  As  in  ref.  [5] 
tire  normal  complete  elliptic  integrals  of  first  and  sec¬ 
ond  kind  are  denoted  K(k)  and  E(k),  respectively. 

Using  eq.  (3)  we  get  the  new  expression  for  the 
hamiltonian 

//  =  8(  1  -u2yll2E(k)lk 

-4*-‘(l  -*2)(1  -u2)'l2K(k).  (5) 

Inserting  eq.  (5)  with  u  =  u(t)  on  the  lhs  of  eq.  (2) 
and  given  by  (3)  on  the  rhs  of  eq.  (2)  we  obtain 

du/dr  =  ±\iry'(\  -  u2)212  -  a'u(  1  -  u 2)  (6) 

with  a'  =  2 ak2E(k)lA,p'  =  20[(2  -  k2)E(k) 


Fig.  1 .  Stationary  one-fluxon  velocity,  um,  versus  length  of 
Josephson  oscillator,  /,  estimated  by  means  of  hamiltonian 
perturbation  theory,  (a)  0  =  0.10,  y  =  0.10,  a  =  0.00,  0.01, 
0.05.  and  0.10,  (b)  a  -  0.02,  y  =  0.05,  0  =  0.00,  0.02, 0.08, 
and  0.20. 


-  2(1  -  k2)K(k)]/A,y'  =  2yk3/A,  and  A  =  2 k*E(k) 
+  (1  -  k2)l2/4K(k).  Here  eq.  (6)  has  the  same  form 
as  the  perturbation  equation  for  the  infinitely  long 
Josephson  transmission  line.  Thus  for  /  -*  00  it  can  be 
shown  that  a’  a, (S'  -»  0,  and  7’  -»  7.  In  the  opposite 
limit,  /  -*■  0,  we  find  a'  -»  5  a,  0'  -*  0,  and  7'  -*  0. 

The  stationary  velocity  u  -  u „  is  obtained  front 
(6)  by  letting  d«/d/  =  0.  The  periodicity  condition 
(4)  is  approximately  valid  for  the  perturbed  sine- 
Gordon  equation  in  the  stationary  state.  Using  (4) 
with  u  =  in  (6)  we  obtain  t<_  as  function  of  /. 

This  function  is  shown  in  fig.  la  and  b  for  different 
values  of  the  parameters  a  and  0.  We  note  that  a  max¬ 
imum  velocity  occurs  for  relatively  small  values  of 
the  ratio  a/0.  In  the  limit  a  =  0  it  can  be  shown  that 
u  -*  1  for  /  -*  0.  For  a  &  0,u  -*  0  for  /  -*■  0  due  to  the 
fact  that  effective  bias,  y',  vanishes  in  this  limit  while 
the  effective  loss,  a',  tends  towards  a  finite  value. 

Fig.  2  shows  the  velocity  as  function  of  the  bias 
for  different  values  of  the  length.  The  curver  inter¬ 
sect  because  the  same  velocity,!/,  occurs  for  two  dif¬ 
ferent  values  of  /  for  certain  values  of  a,  0,  and  7  is 
seen  in  fig.  1.  The  similarity  between  the  curves  and 
the  first  zero  field  step  [4]  in  the  1-V  characteris¬ 
tic  is  due  to  the  fact  that  the  normalized  voltage  V 
=  2rtu^jl  for  /  >  1 .  The  current  1 «  7. 

In  order  to  check  the  validity  of  the  perturbation 
theory  we  have  solved  eq.(l)  numerically  with  the 
static  one-fluxon  solution 

<t>=  2 sin-1  [cn((x  -  l/2)lk,k)]  -  sin-l7 


Fig.  2.  The  bias  current,  7,  versus  the  stationary  one-fluxon 
velocity,  u„,  estimated  by  means  of  hamiltonian  perturba¬ 
tion  theory,  a  =  0.0 1,0  =  0.2,/  =  0.4,  1 .6,  and 
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Fig.  3.  Numerical  solution  of  eq.  (1),  shown  in  terms  of  <px,  with  a  =  0.05,  0  =  0.10,  y  =  0. 10,  and  /  =  1.0.  Inset  shows  the  fluxon 
velocity,  u,  versus  time  r  computed  by  means  of  the  numerical  solution. 


in  the  initial  conditions.  The  method  used  is  a  spec¬ 
tral  (Fourier  transform)  treatment  of  the  space  de¬ 
pendence  together  with  a  leap-frog  scheme  in  time. 
The  resulting  acceleration  of  the  fluxon  is  seen  in 
fig.  3.  The  fluxon  velocity  shown  in  the  inset  con¬ 
verges  towards  the  value  u „  num  =  0.309 1  in  good 
agreement  with  the  predicted  value  M„iPert  =  0.3086 
We  conclude  that  hamiltonian  perturbation  theo¬ 
ry  can  be  used  to  predict  the  stationary  one-fluxon 
velocity  on  the  circular  Josephson  oscillator.  The  re¬ 
sults  will  be  useful  for  the  interpretation  of  experi¬ 
mental  measurements  of  the  l-  V  characteristic  for 
this  device. 


We  are  pleased  to  thank  R.D.  Parmentier,  N.F. 
Pedersen,  A.C.  Scott  and  O.  Skovgaard  for  helpful 
discussions. 
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Standard  methods  of  stochastic  processes  are  used  to  study  the  coupling  of  the  sine-Gordon  sys¬ 
tem  with  a  heat  reservoir.  As  a  result  we  find  thermal  phonons  with  an  average  energy  of  k,T  per 
mode.  The  translational  mode  (zero  mode)  is  found  to  carry  an  average  energy  of  jktT.  This  last 
value  is  just  the  energy  in  the  Brownian  motion  of  a  thermal  soliton.  These  results  are  in  agreement 
with  those  obtained  by  the  use  of  a  statistical-mechanical  description  of  a  dilute  soliton  gas.  Con¬ 
nection  of  the  above  results  with  Josephson  junctions  and  the  extension  of  the  analysis  to  more  gen¬ 
eral  equations  is  also  discussed. 


I.  INTRODUCTION 

The  sine-Gordon  equation  has  recently  been  used  to 
model  several  physical  systems  in  contact  with  a  heat 
reservoir  at  a  given  temperature. 1-9  The  effect  of  the  cou¬ 
pling  of  the  sine-Gordon  system  to  the  heat  reservoir  pro¬ 
vides  a  twofold  mechanism: 

(i)  a  dissipation  of  energy  in  the  system  due  to  an  ener¬ 
gy  flow  from  the  system  to  the  heat  reservoir, 

(ii)  a  disordered  input  of  energy  into  the  system  due  to  a 
flow  back  of  energy  from  the  reservoir. 

A  loss  term  in  the  sine-Gordon  equation  is  then  intrinsi¬ 
cally  connected  to  a  thermal  noise  term,  suggesting  a 
modeling  of  the  interaction  between  the  system  and  the 
reservoir  with  a  driving  stochastic  force  (temperature 
dependent)  in  the  pure  sine-Gordon  equation:1 1-5,7 

•  (I  D 

The  first  term  on  the  right-hand  side  (rhs)  of  Eq.  (1.1)  is 
the  loss  term  representing  the  energy  flow  to  the  reservoir, 
while  the  second  term  is  the  noise  associated  with  a,  giv¬ 
ing  the  disordered  thermal-energy  input  to  the  system. 
The  noise  term  is  assumed  to  be  “white”  both  in  space 
and  time  with  the  autocorrelation  function: 

(n{x,t)n(x',t'))  =  \6a{kBT/E0)8(x  —x')S(t  —t')  .  (1.2) 

Here  (  •  •  •  )  means  ensemble  average,  while  the  constant 
\6a{ktT /E0)  is  determined  by  applying  the  fluctuation 
dissipation  theorem  for  a  soliton  with  small  velocity ',4,7, 10 
iE0  is  the  rest  energy  of  a  soliton  and  is  used  to  fix  the 
scale  of  energy  in  the  system,  kB  is  the  Boltzmann  con¬ 
stant,  and  T  is  the  temperature). 

When  a=0,  n{x,t)=0,  Eq.  (1.1)  reduces  to  the  pure 
sine-Gordon  equation  with  the  exact  soliton  solution 


x  —  vt  I 

exp± 

(l-uJ),/I 

and  Hamiltonian  density 

«  =  -y-[7(^  +*?>+, -«**].  (14) 

Small  oscillations  around  a  ground  state  0O  of  the  system 
are  obtained  by  linearizing  the  sine-Gordon  equation  with 

0  =  0 o+0  >  (1.5) 

this  providing  a  linear  equation  for  0: 

0**— 0h— 0COS0 o=0  >  0  6) 

with  an  associated  energy  density  given  by 

Wph  =  :~(^+^  +  0Wo).  (1.7) 


When  the  ground  state  of  the  system  is  given  by  (1.3),  a 
“zero  mode”  (translational  mode)  is  found  from  Eq.  (1.6). 
In  addition  to  this  mode,  there  exists  a  continuum  set  of 
states  (phonon  modes)  which  satisfy  the  linear  dispersion 
relation:11 

<o2=\+k2  .  (1.8) 

For  practical  applications  to  Josephson  junctions,  it  is  of 
interest  to  include  also  in  the  rhs  of  Eq.  (1.1)  a  constant 
bias  term  tj,  representing  an  ordered  energy  input  into  the 
system  (work  on  the  system).  In  this  case  (in  the  absence 
of  solitons),  phonons  (also  called  “plasmons”)  are  seen 
as  small  oscillations  around  the  ground  state  0() 
=  —  sin-lrj.12 

In  this  paper  we  study  the  effect  of  the  heat  reservoir 
both  on  solitons  and  phonons  by  using  standard  methods 
of  stochastic  processes.  This  will  be  done  in  the  following 
cases. 

In  Sec.  II  we  study  thermal  phonons  in  the  presence  of 
a  static  “exact”  sine-Gordon  soliton.  In  Sec.  Ill  we  in¬ 
clude  a  t]  bias  term  in  the  rhs  of  (1.1)  and  study  the 
thermally  excited  plasmons  around  0O=  —  sin-,Tf  (no  soli¬ 
tons  present  in  the  system). 

In  both  cases  we  find  that  as  long  as  kBT  «£0,  the 
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phonon  modes  have  an  average  energy  of  kgT  per  mode. 
In  the  case  when  a  soliton  is  in  the  system,  however,  we 
And  that  the  corresponding  extra  mode  (zero  mode),  has 
an  average  energy  of  jkBT.  All  these  results  are  in  agree¬ 
ment  with  a  statistical-mechanical  description  of  the 
sine-Gordon  system.  In  Sec.  IV  we  concentrate  on  the  ef¬ 
fect  of  both  the  bias  term  (77)  and  the  heat  reservoir  [a 
and  n(x,t)]  on  the  soliton  motion.  As  result  of  the 
“thermalization"  the  soliton  will  execute  a  Brownian 
motion  with  average  energy  of  \kBT  (the  zero-mode  ener¬ 
gy).  In  Sec.  V  we  relate  the  above  results  to  a  practical 
Josephson  junction  and,  finally,  in  Sec.  VI  we  give  a  short 
summary  of  the  main  results  of  the  paper  including  a 
brief  discussion  of  the  possibility  of  extending  the  analysis 
to  other  equations  of  the  nonlinear  Klein-Gordon  class. 

II.  THERMAL  PHONONS  IN  THE  PRESENCE 
OF  A  SOLITON 

We  consider  as  an  “unthermalized”  system  the  pure 
sine-Gordon  equation 

<!>„-  sin0=O  (2.1) 

and  assume  that  only  a  static  soliton  is  present  (dilute-gas 
limit).  Phonon  modes  are  obtained  as  solutions  of  Eq. 
(1.6)  with  ^0 given  by  (1.3)  with  v—0.  Assuming  i/ik  as 

tM*,<)=/*(x)*'"*',  (2.2) 

we  obtain  from  Eq.  (1.6)  that 

(-3«  +  l— 2sech 2x)/*U)=<4/*(x) .  (2.3) 

As  is  well  known,  (2.3)  admits  a  continuum  set  of  eigen¬ 
functions: 

/*(*)  = - L^eifc,(l+k2)-,/2(k+/tanhx)  (2.4) 

(2  rr)l/2 

together  with  a  zero  mode: 

/»(*)  = -^sech*  (2.5) 

which  restores  the  translational  symmetry  broken  by  the 
introduction  of  the  soliton  into  the  system  (Goldstone 
mode).11  Equations  (2.4)  and  (2.5)  together  form  a  com¬ 
plete  set  of  orthonormal  eigenfunctions: 

/_ ^  f£(x)dx  =  1,  /_+^  fb(x)fk{x)dx  =0 
f*~fl(x)fk\x)dx  =6(* -k')  ,  (2.6) 

/»(x)/»(x')  +  J+yUx)fk{x')dk=b(x  -x') 

where  *  in  the  superscript  means  complex  conjugate. 

By  coupling  the  sine-Gordon  system  with  the  heat 
reservoir,  we  change  Eq.  (2.1)  into  Eq.  (1.1).  Thermal 
phonons  are  then  found  to  satisfy 

tl>XJ-}f>„-ifrcos<fio=aipl+nlx,t)  (2.7) 

for  which  the  general  solution  can  be  expa.  Jed  in  terms 
of  the  complete  set  (2.6)  as 


*klx,t)=  2  AkU)fk(x)  +  Ab(t)fb(x)  (2.8) 

*  -0 

(here  we  assume  the  system  to  be  in  a  box  of  length  L, 
and  then  let  L— ►  »).  Substituting  (2.8)  in  (2.7)  and  using 
Eq.  (2.3),  we  obtain 

2  lAk.„fklx)+aAklfk{x)+Ail<0lfk{x )J 
»-o 

+  (Abn+aAbl)fb(x)=  —n(x,t)  .  (2.9) 

Equation  (2.9)  is  easily  studied  once  projected,  respective¬ 
ly,  along  the  fklxYs  and  the  fb(x)  eigenfunctions,  this 
giving  [using  (2.6)] 

Ak,u+aAkl  +  Akc)k=(k{t )  (2.10) 

and 

Ab,tt+aAb'i=(k(t)  (2.11) 

with 

eb(t)=  —  f  fh(x)n(x,t)dx  , 

"  —  to 

ek(t)=  —  f  fl{x)n(x,t)dx  . 

"  —  m 

By  using  (1.2)  and  (2.6)  we  find,  for  the  autocorrelation 
function  RtU  — /')  and  for  the  power  spectrum  S,{a)  of 
the  normal  processes  ek(t)  and  eb(t),  that 

R, k(t-f)~RuU~t’) 

=  16 a{kBT/E0)bU-i'),  (2.12) 

S, h(a)=Su(o})=l6a(kBT/E0) .  (2.13) 

Equations  (2.10)  and  (2.11)  are  then  integrated  by  the 
standard  theory  of  stochastic  processes,10  giving  the  fol¬ 
lowing  expressions  for  the  power  spectrum  of  Ak(t)  and 
Ak.Mh 

SAA<o)=l6a{kBT/E0)—1 - ^ - j-j-  .  (2.14) 

*  (or— <w*r+cr« 

S44i(<u)  =  w2S^(<u)  .  (2.15) 

[S^(tu)  and  S^lw)  are  obtained  from  (2.14)  and  (2.15) 
with  the  substitution  uk  =0.]  If  we  assume  ergodicity,  the 
time  averages  of  the  processes  |  i4tU)|2  and  M*,,(f)|2 
are  evaluated  as 

<M*(0|2>=r,4(0) 

=  jy  j^SAt{a>)  =  UkB T/o>l E0 ) ,  (2.1 6) 

(|^.,(f)|2>=^t|(0) 

r  +  («)  =  8<*,7Y£0) ,  (2.17) 

*  “  •  2rr 

where  contour  integration  has  been  used  in  evaluating  the 
integrals  in  (2.16)  and  (2.17).  In  the  same  way  we  obtain 
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for  < 

<M»,,('>|J>=«*«7Y£o).  <2. 18) 

From  (1.7)  and  (2.8)  we  have  that  the  average  energies  of 
the  Jkth  phonon  mode  Ak{t)fk(x)  and  of  the  translational 
mode  Ab(t)fk(x)  are  given,  respectively,  by 

{Hk)=l*[{  M*.»(,)|2>+^<  M*{')|2>]  <2.19) 

and 

(2.20) 

Using  (2.16)— (2.18),  we  finally  have 

(Hk)  —kBT  (2.21) 

and 

(Hk)  =  \kaT  (2.22) 

in  complete  agreement  with  the  classical  statistical- 

mechanics  analysis  of  a  dilute  soliton  gas  derived  in  Ref. 

8. 

III.  THERMAL  PHONONS  IN  THE  ABSENCE 
OF  SOLITONS  WITH  BIAS 

In  this  section  we  consider  the  unthermalized  system  to 
be  the  pure  sine-Gordon  system  of  finite  length  L,  with  a 
constant  driving  force  77  <  1: 

sin d=i?  . 


with 


and  (2/L),/J  being  just  a  normalization  factor  (for  0 
it  should  be  read  as  L~in).  Substituting  (3.7)  in  (3.S)  and 
applying  to  both  sides  of  the  equation  the  pro¬ 
jection  operator  f^cos(k„x)dx,  we  obtain 

A*«+aAm'l  +  Ul-V1)i/2+k;]Am=emU) .  (3.9) 

where 

e,(t)=(2/L)l/J  fQ  n{x,t)cos(k,,x)dx  .  (3.10) 

Using  (1.2),  we  find  for  the  autocorrelation  function 
and  the  power  spectrum  S,^(t o)  the  same  ex¬ 
pression  as  in  (2.12)  and  (2.13).  By  identifying 

[  ( 1  —  T)1 ) 1  n + kl  ]  with  a>l ,  we  see  that  Eq.  (3.9)  in  the  lim¬ 
it  L-*co  coincides  with  Eq.  (2.0),  and  therefore,  follow¬ 
ing  the  same  analysis  of  the  preceding  section,  we  obtain 
that  the  average  energy  per  phonon  mode  is 

( H„)=kBT .  (3.11) 

No  zero-mode  energy  is  present  in  this  case,  due  to  the  ab¬ 
sence  of  the  soliton  in  the  system.  Finally,  we  remark 
that  the  above  results  do  not  depend  on  the  particular 
boundary  condition  (3.3)  used  (we  could  have  used  generic 
periodic  boundary  conditions)  as  well  as  on  smallness  re¬ 
quirements  of  a  and  rj.  The  only  approximation  that  has 
been  made  in  obtaining  (2.21),  (2.22),  and  (3.11)  is  vhe 
linearization  procedure,  which  is  justified  if 

(3.12) 


(3.1)  kBT «E0  , 


Phonon  modes  are  seen  as  small  oscillations  around 
the  classical  ground  state 

<fi0=  —  sin-,7j  ,  (3.2) 

satisfying  the  boundary  conditions 

0«.*(O,f)=^,.,a,f)=O  (3.3) 

(no  solitons  are  present  in  the  system).  The  thermalized 

system  is  obtained  from  Eq.  (3.1)  by  adding  to  the  rhs  the 
term  a#,  +  n(x,t)  with  n(x,t)  given  as  in  (1.2): 

dju-^n—  sin^=i/+a^,+nU,f) .  (3.4) 

Thermal  phonons  are  then  solutions  of  the  following  sto¬ 
chastic  equation: 

■■  (3.5) 


as  appears  evident  from  Eqs.  (2.16)  and  (2.17). 

IV.  BROWNIAN  MOTION  AND  DIFFUSION 
CONSTANT  OF  A  THERMAL  SOLITON 

We  now  concentrate  on  the  effect  of  the  a,  r\,  and 
n(x,f)  terms  in  Eq.  (3.4)  on  the  soliton  motion  (here  a  sol¬ 
iton  is  a  2?r-kink  jump  from  —  sin-1?/  to  2?r— sin-1?/). 
We  assume  17 /a  and  kBT/E0  to  be  small.  By  introducing 
the  momentum 

f>=  — 7  f  +  "°tx$idx  ,  (4.1) 

and  differentiating  with  respect  to  time,  we  obtain 

^■  =  -aP+iin)+eU) ,  (4.2) 


When  a— 0,  n(x,l)= 0,  these  phonons  are  just  classical 
Klein-Gordon  modes  with  energy  given  by 

+*?+#*( l-i72)1/2J  •  <3.6) 

The  general  solution  of  Eq.  (3.3)  satisfying  the  boundary 
conditions  (3.3)  is  of  the  form 

^=(2/L),/,2^U)cos(M)  (3.7) 


where  we  have  used  Eq.  (3.4)  to  eliminate  the  term  and 
have  defined  e(r)  as 

c( r)  =  j J"+  $x(x,t)n(x,t)dx  .  (4.3) 


Neglecting  the  noise  term,  Eq.  (4.2)  describes  the  “power 
balance"  motion  of  a  2ir-kink  with  velocity11 


«o  =  ±  1  + 


(4.4) 
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and  momentum 


Po= 


up 

(1  -u\)'n 


(4.5) 


The  noise  term  n  ( x,t )  in  Eq.  (4.2)  introduces  fluctuations 
in  the  momentum  and,  from  (4.5),  in  the  velocity  of  the 
kink.  Such  fluctuations  are  readily  evaluated  by  standard 
techniques  (for  details  see  Ref.  7),  giving  the  following  for 
the  power  spectrum  of  the  process  Au  (/)  =  u  (r)  —  u0: 


(, -ul)in 

SBJa)^2a(kBT/E0) - - - r—  .  (4.6) 

co2+a2 

By  assuming  u0  «  1,  Eq.  (4.6)  reduces  to 

Su(<o)  =  2a(kBT/E0)-~-^  (4.7) 

io2+a 2 

from  which  we  obtain 

<w2>  =  f+’‘^LSu(co)  =  (kBT/EQ)  .  (4.8) 

The  time  average  of  the  kinetic  energy  in  the  Brownian 
motion  of  the  soliton  is  then  given  by 

<£•„,)  =  ±E0(u2)  =  ±kBT  (4.9) 

as  expected  from  soliton  statistical-mechanics  theory.* 
Finally,  from  Eq.  (4.6)  a  diffusion  constant  D  for  the 
2rr-kink  motion  is  derived: 

D  =  ~(kBT/a)  ,  (4.10) 


which  is  just  the  usual  Einstein  diffusion  constant  for  the 
Brownian  motion  of  a  particle  in  a  viscous  medium  (this 
is  a  further  confirmation  of  the  particlelike  nature  of  the 
soliton).  (See  also  Ref.  3.) 


V.  THERMAL  SOLITONS  AND  PHONONS 
IN  JOSEPHSON  JUNCTIONS 

We  will  now  relate  the  foregoing  sections  to  a  real  de¬ 
vice  as  the  Josephson  junction.  We  will  find  the  orders  of 
magnitude  of  the  quantities  of  interest  and  see  if  the  as¬ 
sumption  made  in  the  above  analysis  holds  for  Josephson 
junctions. 

The  fluxon-rest  energy  (in  laboratory  units)  for  a 
Josephson  junction  is 

E0  =  Au0 = 8M. jJL  /( 2e )  ,  (5.1) 

where  J  is  the  maximum  Josephson  current  density,  L  is 
the  length  of  the  junction,  and  e  is  the  electron  charge,  kj 
in  (5.1)  is  the  Josephson  penetration  depth  given  by 

kJ=(fi/2eti0dJ)ul  (5.2) 

where  d  is  the  magnetic  thickness  of  the  oxide  layer 
(2kL  +  /<>)>  and  ji0  the  vacuum  permeability.  From  (5.1) 
we  have  that  for  a  typical  Josephson  junction 

(*,7y£0)«10-4-10-\  (5.3) 

which  justifies  the  assumption  (3.12)  made  in  the  analysis. 


The  Josephson  plasma  frequency  is 

6>/.  =  (2e/rQ/e0e,fl)l/J  ,  (5.4) 


where  er  and  e0  are,  respectively,  the  relative  dielectric 
constant  of  the  oxide  layer  and  the  dielectric  constant  of 
the  vacuum,  while  »0  is  the  thickness  of  the  oxide  layer. 
For  a  plasmon  described  in  Sec.  Ill  we  have  that  the  split¬ 
ting  of  the  energy  level  is 


£pi  =**>/.*>. 


(5.5) 


with  a),,  given  by  [(1—  7)2)l/2  +  {nir/L)2],/2.  We  have 
then  that  the  ratio  kBT /£*j  is  of  the  order  of  magnitude 
1  —  10,  i.e.,  the  quantum  energy  levels  are  separated  by  a 
quantity  comparable  with  kBT.  To  have  a  rough  estimate 
of  the  energy-level  separation  for  a  fluxon,  we  can  use  the 
analogy  of  a  particle  in  a  box.  This  gives 


i rffi  n 
2 M0  L2 


\0~**op  , 


(5.6) 


.e.,  for  a  fluxon  the  separation  in  the  energy  levels  is 
smaller  than  kBT  by  a  factor  of  the  order  10~J— 10-4. 
This  numerical  manipulation  indicates  that  for  a  typical 
Josephson  junction  fluxon  quantitation  is  not  necessary, 
while  it  is  necessary  for  plasmons  (£J,  being  of  the  same 
order  of  magnitude  as  kBT).  In  Ref.  6  the  effects  of 
quantum  plasmons  on  the  fluxon  motion  have  been  calcu¬ 
lated.  It  turns  out  that  they  are  several  orders  of  magni¬ 
tude  smaller  than  the  direct  influence  of  the  thermal 
reservoir  on  the  soliton  evaluated  in  this  paper,  and  there¬ 
fore,  in  our  context,  completely  negligible. 

Finally,  in  closing  this  section  it  is  worth  noting  that  if 
kBT/E0  is  very  small,  a  statistical-mechanical  description 
of  fluxons  in  Josephson  junctions  is  meaningless.  Howev¬ 
er,  the  method  used  in  the  preceding  section  is  still  useful 
to  study  the  interactions  between  plasmons  and  fluxons. 
(See  Ref.  7  for  the  case  of  Josephson  oscillators.) 


VI.  CONCLUSION 

It  has  been  shown  that  the  effect  of  a  thermal  reservoir 
on  the  sine-Gordon  system  can  be  studied  by  using  stan¬ 
dard  methods  of  stochastic  processes.  Both  phonons  and 
solitons  are  found  to  be  thermalized  in  a  way  that  the 
phonons  will  have  an  average  energy  of  kBT  per  mode, 
while  solitons  will  have  an  energy  of  \kBT.  These  results 
are  in  agreement  with  those  obtained  by  using  a 
statistical-mechanics  approach  for  a  "dilute”  solution 
gas.8  The  main  assumption  used  in  our  derivation  has 
been  kBT «E0  (to  justify  the  linearization  procedure). 
Second-order  effects  [in  the  small  quantity  (kBT/E0)\, 
such  as  interaction  between  phonon  modes  and  solitons,4-4 
have  been  neglected  therefore.  Finally,  in  closing  this  pa¬ 
per  we  wish  to  point  out  that  in  spite  of  the  particularity 
of  the  model  used,  the  results  obtained  are  sufficiently 
general  to  be  extended  to  other  equations  of  the  nonlinear 
Klein-Gordon  class,  such  as  ^4,  double  sine-Gordon,  etc. 
As  a  matter  of  fact,  the  only  difference  in  the  analysis  will 
be  the  presence  of  additional  bound  states  in  the  linear 
phonon  eigenvalue  problem.  By  following  arguments 
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similar  to  those  used  in  the  present  paper,  however,  it  is 
easily  shown  that  each  of  these  additional  bound  states 
carries  a  thermal-average  energy  of  kBT.  This  energy  will 
not  increase  the  energy  of  the  center  of  mass  of  a  soliton- 
like  solition  of  these  more  general  models,  but  it  will  in¬ 
crease  the  energy  of  the  internal  degrees  of  freedom 
(motion  around  the  center  of  mass)  of  these  excitations.14 
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Experimental  measurements  of  current-voltage  structure  and  emitted  Y-band  radiation  in  applied 
magnetic  field  from  overlap-geometry  Josephson  tunnel  junctions  of  normalized  length  about  2  are 
compared  with  numerical  simulations  obtained  with  the  use  of  a  perturbed  sine-Gordon  model.  The 
simulations  furnish  the  current  and  field  dependence  of  the  oscillation  configuration,  from  which 
can  be  calculated  average  voltages,  frequencies,  and  power  spectra.  Simulation  and  experimental  re¬ 
sults  are  in  good  agreement  with  regard  to  the  lobe  structure  of  the  height  of  the  first  zero-field  step 
and/or  second  Fiske  step  in  magnetic  field  and  the  field  dependence  of  the  radiation  frequency 
within  the  various  lobes,  including  details  such  as  hysteresis  between  lobes.  The  simulations  predict 
an  alternation  of  the  dominant  frequency  component  with  increasing  field  that  accounts  well  for  the 
experimental  observations.  The  usefulness  and  limitations  of  cavity-mode  analyses,  both  single¬ 
mode  and  multimode,  are  evidenced  by  comparison  with  the  simulation  results. 


I.  INTRODUCTION 

Fluxon  (soli ton)  propagation1  has  been  by  now  estab¬ 
lished  as  the  essential  physical  mechanism  underlying 
many  of  the  observed  experimental  properties  of  long 
Josephson  tunnel  junctions.  In  particular,  the  appearance 
of  both  zero-field  steps  (ZFS’s)  and  Fiske  steps  (FS*s)  in 
the  current-voltage  ( /-  V)  characteristics  of  such  junctions 
and  the  emission  of  microwave  radiation  from  junctions 
when  biased  on  these  steps  seem  to  be  explainable  in  terms 
of  fluxon  dynamics.2-5  A  number  of  different  ap¬ 
proaches  have  been  employed  in  the  literature  to  account 
for  the  available  experimental  observations.  These  include 
perturbative  expansions  of  the  basic  soliton  equation  in¬ 
volved  (sine-Gordon  equation),6  analytic  extensions  of 
small-junction  theory  (cavity-mode— interaction  analy¬ 
ses),7  and  mechanical  analog*  and  digital  computer9  simu¬ 
lations  (the  references  cited  are  intended  to  be  representa¬ 
tive,  not  exhaustive).  Moreover,  direct  dynamic  measure¬ 
ments  at  the  single-fluxon  level  are  beginning  to  appear  in 
the  literature.10-" 

The  perturbative  approach  is  most  suited  for  studying 
the  behavior  of  low-order  steps  on  ary  long  junctions, 
inasmuch  as  the  usual  point  of  departure  here  consists  of 
the  exact  analytic  solutions  of  the  sine-Gordon  equation 
on  the  infinite  spatial  interval.  Multimode  extensions  of 
small-junction  theory,  on  the  other  hand,  should  presum¬ 
ably  be  most  appropriate  for  relatively  short  junctions. 
For  junctions  which  are  neither  very  long  nor  very  short, 
one  would  expect  a  priori  that  neither  of  these  two  ap¬ 
proaches  could  be  counted  on  to  give  reliable  results.  In 
such  cases,  direct  simulation  would  seem  to  be  indispens¬ 
able. 

The  present  paper  is  an  attempt  to  elucidate  further  and 


in  more  detail  just  this  case,  viz.,  the  dynamics  underlying 
the  behavior  of  intermediate-length  Josephson  junctions. 
To  this  end  we  compare  the  results  of  new  experimental 
measurements  of  I-  V  structure  and  microwave  emission 
in  magnetic  field  with  the  results  of  detailed  numerical 
simulations.  The  agreement  that  emerges  is  quite  con¬ 
vincing.  For  simplicity,  attention  is  focused  primarily  on 
the  first  zero-field  step  (ZFSl)  and  on  the  second  Fiske 
Step  (FS2)  in  junctions  of  normalized  length  of  about  2. 
Since  the  voltage  positions  of  ZFSl  and  FS2  approximate¬ 
ly  coincide,  we  refer  in  the  following  to  the  step 
ZFS1/FS2.  The  numerical  simulations  are  compared  also 
with  approximate  analytic  results,  and  the  usefulness  and 
limitations  of  the  latter  are  clarified. 

The  paper  is  structured  as  follows:  Section  II  contains 
a  description  of  the  mathematical  model  used  and  the 
techniques  employed  in  its  analysis.  The  results  of  this 
analysis  are  presented  in  Sec.  III.  The  experimental  mea¬ 
surements  performed  are  described  and  discussed  in  Sec. 
IV.  Finally,  Sec.  V  contains  our  concluding  comments. 

II.  MATHEMATICAL  MODEL 
AND  COMPUTATION  TECHNIQUES 

The  mathematical  model  studied  is  the  perturbed  sine- 
Gordon  equation,9 

4>x*  —  —  sin<£=a<£,  —  —y  ,  (la) 

**«y)=<fc,a, »>  =  ■*)  .  Ub) 

appropriate  to  an  overlap-geometry  junction."  Here,  4>  is 
the  usual  Josephson-phase  variable,  x  is  distance  along  the 
junction  normalized  to  the  Josephson  penetration  depth 
Kj,  l  is  time  normalized  to  the  inverse  of  the  Josephson 
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plasma  (angular)  frequency  <o0,  and  subscripts  denote  par¬ 
tial  derivatives  (see  Ref.  9  for  details  of  the  normaliza¬ 
tions).  The  model  contains  five  parameters:  a,  0,  y,  L, 
and  17.  The  term  in  a  represents  shunt  (quasiparticle)  loss. 
The  0  term  models  dissipation  due  to  the  surface  resis¬ 
tance  of  the  superconducting  films.  The  constant  y  mea¬ 
sures  a  uniformly  distributed  bias  current  normalized  to 
the  maximum  zero-voltage  Josephson  current.  The  con¬ 
stant  17  is  a  normalized  measure  of  the  external  magnetic 
field  which  determines  the  boundary  conditions  [Eq.  (lb)] 
at  the  two  ends  of  the  junction  of  normalized  length  L. 
In  this  numerical  study  the  dissipative  and  length  parame¬ 
ters  were  fixed  at  a =0.05,  0=0.02,  and  L  =  2.  These 
were  chosen  to  be  representative  of  typical  physical  values 
without,  however,  modeling  any  one  specific  junction. 
The  bias  current  and  magnetic  field  parameters  were 
varied  in  the  ranges  0<y  <  1,  0<ij  <6. 

Equations  (1)  were  integrated  from  initial  conditions  ei¬ 
ther  similar  to  Eqs.  (2)  of  Ref.  13  or  (more  often)  using 
the  final  p  and  p,  distributions  from  runs  at  nearby  points 
in  parameter  space.  The  integration  was  carried  out  using 
the  implicit  finite-difference  method  described  in  detail  in 
Ref.  9,  with  space  and  time  intervals  set  to  0.02  and  0.01, 
respectively.  Numerical  accuracy  and  stability  were  veri¬ 
fied  by  halving  and  doubling  these  intervals.  During  each 
run  the  time-averaged  value  of  the  voltage  at  the  two  ends 
of  the  junction  and  the  power  spectrum  of  the  voltage  at 
the  left  (x=0)  end  were  calculated.  Here,  voltage  is  de¬ 
fined  as  p„  which  represents  the  physical  voltage  normal¬ 
ized  to  fko0/2e,  where  ft  is  Planck's  constant  divided  by 
2ir,  and  e  is  the  magnitude  of  the  electronic  charge. 
These  quantities  were  calculated  over  an  integral  number 
of  oscillation  periods  during  the  last  approximately  50 
normalized  time  units  of  each  run.  The  power  spectra 
were  calculated  by  means  of  a  fast  Fourier  transform  us¬ 
ing  a  simple  rectangular  window.14  The  values  of  (p,) 
were  calculated  both  from  the  elementary  definition  of 
average  and  the  zero-frequency  components  of  the  power 
spectra. 

Two  checks  were  employed  to  assure  that  the  average 
voltages  and  the  power  spectra  were  calculated  over 
steady-state,  not  transient,  dynamic  configurations:  (i) 
The  values  of  (p, )  at  the  two  junction  ends  were  com¬ 
pared.  Physically  and  mathematically,  the  time-averaged 
voltage  must  be  constant  along  the  length  of  the  junction, 
(ii)  The  quantity  ( P,)/2irf ,  where  /  is  the  fundamental 
oscillation  frequency,  was  calculated.  From  the  Joseph¬ 
son  frequency  relation,  this  quantity,  in  steady  state,  must 
be  an  integer  whose  value  (1  or  2  in  the  present  case)  de¬ 
pends  upon  the  type  of  oscillation  present.*  If  either  of 
these  conditions  was  not  satisfied  to  within  specified  lim¬ 
its  the  time  duration  of  the  run  was  increased. 

III.  NUMERICAL  RESULTS 

Figure  1  shows  the  magnetic  field  diffraction  pattern  of 
ZFS1/FS2.  In  this  figure,  circles  and  diamonds  are  the 
numerically  computed  top  of  the  step.  Circles  were  calcu¬ 
lated  by  increasing  y  at  constant  17;  diamonds  were  deter¬ 
mined  by  varying  i\  at  constant  y  (as  will  be  seen  later,  the 
difference  is  significant).  Small  arrows  near  the  diamonds 
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FIG.  1.  Magnetic  field  diffraction  pattern  of  ZFSI/FS2. 
Circles:  step  top  calculated  numerically  at  constant  tj.  Squares: 
step  bottom  calculated  numerically  at  constant  r).  Diamonds: 
step  top  calculated  numerically  at  constant  y.  Arrows  near  dia¬ 
monds  indicate  direction  of  field  variation.  Solid  curve  (marked 
1  and  2):  Enpuku  el  at.  multimode  theory.  Dashed  curve:  Ku- 
lik  single-mode  theory.  Enpuku  and  Kulik  curves  coincide  at 
the  bottom  of  the  first  lobe.  Insets  show  approximate  dynamic 
trajectories  in  the  various  tj  regions:  solid  lines  are  fluxons  or 
antifluxons;  dashed  lines  are  plasma  waves. 


indicate  the  direction  in  which  17  was  varied  near  the  max¬ 
imum  point.  Beyond  the  maximum  points,  the  junction 
switched  to  a  different  dynamic  state,  most  often  to  the 
McCumber-Stewart  hysteresis  curve  at  the  corresponding 
value  of  y.  This  fact  was  evidenced  by  abrupt  changes  in 
the  value  of  (p,),  the  voltage  waveform,  and  the  corre¬ 
sponding  power  spectrum.  In  zero  magnetic  field,  the 
form  of  the  McCumber-Stewart  hysteresis  curve  can  be 
approximated  by15,2 

y=4aE{k)/irk  ,  (2a) 

(p,)  =ir/kK{k) ,  (2b) 

where  K(fc)  and  E(k)  are,  respectively,  the  complete  el¬ 
liptic  integrals  of  the  first  and  second  kinds  of  modulus  k. 
From  Eqs.  (2)  it  follows  that  for  k—*0,  y-*a{p,),  i.e., 
the  McCumber-Stewart  curve  approaches  asymptotically 
the  Ohmic  line,  whereas  for  k-*  1,  {p,  >— -0,  and 
y—*4a/ir.  Equations  (2)  continue  to  hold  as  a  rough  ap¬ 
proximation  even  in  the  presence  of  magnetic  field,  at 
least  for  rj  <  1 . 

The  squares  in  Fig.  I  are  the  numerically  computed 
bottom  of  the  step,  calculated  by  decreasing  y  at  constant 
77.  As  will  become  clear  from  Figs.  2—6,  it  is  possible  to 
determine  numerically  the  precise  bottom  of  the  step  in 
the  first  lobe  of  the  diffraction  pattern,  but  this  is  no 
longer  the  case  in  the  second  lobe.  Beyond  the  minimum 
points,  the  solution  followed  the  McCumber-Stewart 
curve  down  for  a  certain  distance,  after  which  it  switched 
abruptly  to  a  static,  zero-voltage  state.  The  lowest  y  value 
for  which  such  switching  was  observed  was  0.06 
^y<0.07,  which  is  consistent  with  the  value 
4a  /ir*  0.064  estimated  from  Eq.  (2a). 

The  dashed  curve  in  Fig.  1  has  been  calculated  from  the 
Kulik  theory14  for  FS2.  The  input  parameter  for  this 


2642 


M.  P.  SOERENSEN  et  al. 


30 


theory  is  the  quantity  Z„=L1QK/n2n2,  where  n  is  the 
step  number  (2,  in  our  case)  and  Q„  is  the  quality  factor 
of  the  nth  mode.  Following  Enpuku  et  al.  (Ref.  7)  Q„  is 
defined  by 

_L._iSL+ JLE§.  (3) 

Q,  nir  L 

Insertion  of  parameter  values  thus  yields  Z2  =  1.2867. 
The  Kutik  theory  gives  the  maximum  height  of  the  step 
above  the  Ohmic  current.  Accordingly,  to  compare  with 
numerical  (or  experimental)  results  it  is  necessary  to  add 
this  component  to  the  Kulik  value.  Since,  from  our  nu¬ 
merical  results,  the  top  of  FS2  is  at  ($, )  =3.1,  a  constant 
(independent  of  17)  value  of  a{$, )  =0.155  has  been  added 
to  the  Kulik  component  in  drawing  the  dashed  curve  in 
Fig.  1. 

The  Kulik  theory  is  seen  to  predict  the  maximum 
points  of  the  second  lobe  up  to  17 «  4.  However,  for  the 
maximum  points  of  the  first  lobe  and  for  the  second  lobe 
for  17  >  4,  the  theory  fails.  To  predict  these  results  more 
modes  must  be  included  in  the  computations.  Following 
Enpuku  et  al.1  we  have  used  five  modes  with  the  quality 
factor  Q„  given  by  Eq.  (3)  with  n=  1,2 .....  5.  We  have 
solved  Eqs.  (10)— (12)  in  Ref.  7  by  means  of  a  standard 
routine.17  The  results  are  shown  as  the  solid  curve  in  Fig. 
1.  The  Enpuku  theory  is  seen  to  predict  the  maximum 
points  of  the  first  lobe  very  well.  The  maximum  values 
on  the  second  lobe  up  to  77  =  3.3  are  also  in  agreement. 
However,  between  i|  —  3.1  and  5.5  the  Enpuku  theory 
predicts  two  curves  for  the  maximum  values.  The  curves 
are  marked  1  and  2  in  accordance  with  a  major  contribu¬ 
tion  to  the  solution  from  the  first  and  second  cavity 
modes  respectively.  Curve  2  lies  close  to  the  Kulik  curve, 
in  agreement  with  the  fact  that  this  latter  curve  was  com¬ 
puted  by  means  of  mode  2  exclusively.  The  numerically 
computed  maximum  values  (circles)  agree  with  the  upper 
curves  (i.e.,  curve  2  in  the  interval  17  =  3. 1—4.4  and  curve 
1  in  the  interval  17=4.4—5.5).  Just  above  the  lower  curves 
(i.e.,  curve  1  for  17  =  3.1—4.4  and  curve  2  for  ij=4.4— 5.5) 
the  computer  results  exhibit  changes  in  the  contents  of 
cavity  modes  from  mode  l  to  mode  2  above  curve  1  and 
vice  versa  above  curve  2.  The  accompanying  hysteresis 
phenomena  are  discussed  below. 

Figures  2—6  depict  five  vertical  (constant-17)  slices 
through  the  diffraction  pattern.  Figure  2(a)  shows  the 
current  dependence  of  the  power  levels  of  the  first  two 
Fourier  harmonics  of  the  voltage  at  the  left  end  of  the 
junction  in  zero  magnetic  field.  Power  levels  are  given  by 
10  In  |  A  | 2  -1-  100,  A  being  the  voltage  Fourier  com¬ 
ponent.  The  ac  components  are  thus  arbitrarily  normal¬ 
ized  to  10-,°.  Since  no  loading  effects  are  included 
in  the  model,  all  power  levels  calculated  should  be  con¬ 
sidered  ideal,  available  values.  Figure  2(b)  shows  the 
current  dependence  of  the  first  harmonic  frequency  f\  of 
the  oscillation.  Since  average  voltage  and  frequency  are 
proportional  (in  steady  state)  through  the  relations 
)=4ir/1=2ir/1,  Fig.  2(b)  is  effectively  the  current- 
voltage  characteristic  of  the  step.  We  have  chosen  to  plot 
this  characteristic  in  frequency  rather  than  in  voltage  be¬ 
cause  in  the  laboratory,  frequency  can  be  measured  much 
more  precisely  (although  perhaps  less  easily)  than  voltage. 
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FIG.  2.  Section  through  diffraction  pattern  at  17=0.  (a) 
Power  levels  of  first  two  Fourier  voltage  components  st  x  —  0, 
normalised  to  I.Ox  I0"'°.  (b>  First-harmonic  frequency. 
Voltage  waveform  at  x=0  and  corresponding  power  spectrum 
are  shown  for  (c)  y =0.60,  (d)  y  =  0.30,  (e)  y  =0. 15,  and  (fi 
y  =0.13. 
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Figures  2(c) — (f)  show  the  voltage  waveform  and  part  of 
the  corresponding  power  spectrum  at  the  four  points  indi¬ 
cated  by  arrows  in  Fig.  2(b).  A  comparison  of  Figs. 
2(a)— 2(e)  with  results  from  the  literature  [see,  in  particu¬ 
lar,  Fig.  17  of  Ref.  9  and  Fig.  3(a)  of  Ref.  2]  leaves  no 
doubt  that  the  oscillation  depicted  in  Fig.  2  is  the  fluxon 
oscillation  corresponding  to  ZFS1.  Finally,  a  comparison 
of  Figs.  2(e)  and  2(0  makes  clear  why  it  is  possible  to 
determine  precisely  the  bottom  of  the  step  in  the  first  lobe 
of  the  diffraction  pattern:  At  a  certain  value  of  the  bias 
current  (0.14 < y  <0.15  in  the  present  case)  the  /,  com¬ 
ponent  of  the  oscillation  abruptly  disappears,  and  a  new 
oscillation  evolves  for  which  the  dominant  component  is 
at/2. 

Figure  3  shows  a  similar  vertical  slice  through  the  dif¬ 
fraction  pattern  at  a  point  near  the  right-hand  extremity 


CURRENT  Y 
«D 


CURRENT  Y 
(W 


Time  t  Frequency  f 


(c) 


(•I 

FIG.  3.  Section  mrough  diffraction  pattern  at  77  =  1.5.  (a) 
Power  levels  of  first  two  Fourier  voltage  components  at  x=0, 
normalized  to  ^=1.0xl0-l°.  (b)  First-harmonic  frequency. 
Voltage  waveform  at  jc=0  and  corresponding  power  spectrum 
are  shown  for  (c)  y =0.34,  (d)  y =0.25,  and  (e)  y  =0. 19. 
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FIG.  4.  Section  through  diffraction  pattern  at  (a) 

Power  levels  of  second  Fourier  voltage  component  at  x =0,  nor¬ 
malized  to  ^?=  1.0X  I0~'°  (P,  is  absent),  (b)  Second-harmonic 
frequency.  Voltage  waveform  at  x=0  and  corresponding  power 
spectrum  are  shown  for  (c)  y  =0.30  and  (d)  y  =0.16. 
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of  the  Tint  lobe  (ij  *■  1.3).  Here,  the  behavior  of  the  oscil¬ 
lation  near  the  top  of  the  step  is  qualitatively  quite  similar 
to  that  shown  in  Fig.  2.  As  y  is  decreased,  however,  a 
marked  difference  is  seen:  The  ratio  of  the  / j  power  Pj 
to  the  fi  power  Pi  steadily  decreases  until,  near  the  bot¬ 
tom  of  the  step,  P2  in  fact  dominates.  Since  /2  is  the 
dominant  frequency  of  FS2,  the  nature  of  the  step  at 
1.3  may  be  described  as  ZFSl-like  near  the  top  and 
FS2-like  near  the  bottom.  It  should  be  noted,  however, 
that  the  transition  between  the  two  oscillation  configura¬ 
tions  is  here  relatively  gradual,  rather  than  abrupt. 

Figure  4  shows  a  vertical  slice  at  17  =  2.5,  in  the  left- 
hand  side  of  the  second  lobe  of  the  diffraction  pattern. 
Here  the  oscillation  is  purely  FS2-like  over  the  entire  ex¬ 
tent  of  the  step.  The  power  Pi  is  lost  in  the  background 
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FIG.  3.  Section  through  diffraction  pattern  at  r/  =  3A.  (a) 
Power  levels  of  second  Fourier  voltage  component  at  x=0,  nor¬ 
malized  to  4? “ I.OX  10“ 10  (P,  it  absent  or  very  small),  (b) 
Second-harmonic  frequency.  Voltage  waveform  at  x=Q  and 
corresponding  power  spectrum  are  shown  for  (c)  y  =  0.42  and  (d) 
r=O.I6. 


noise  (the  noise  here  is  of  numerical,  not  physical,  origin), 
and  P2  is  the  dominant  component.  The  fact  that  the 
next-highest  component,  at  /  =  2/2,  lies  20  dB  or  more 
below  P2,  explains  why  the  Kulik  theory,  which  assumes 
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FIG.  6.  Section  through  diffraction  pattern  at  ij«=  4.6.  (a) 
Power  levels  of  first  two  Fourier  voltage  components  at  *=0, 
normalized  to  4?"  l.Ox  IQ'10,  (b)  Second-harmonic  frequency. 
Voltage  waveform  at  x  =  0  and  corresponding  power  spectrum 
are  shown  for  (c)  y =0.32,  (d)  y=0.28,  and  (e)  y  =0.14. 
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a  single  mode  at  f  =f2,  gives  a  reasonable  description  of 
the  oscillation  in  this  region,  as  seen  in  Fig.  1.  The  diffi¬ 
culty  in  determining  precisely  the  bottom  of  the  step  in 
the  second  lobe  is  now  also  apparent:  f2  is  the  dominant 
component  of  both  FS2  and  the  McCumber-Stewart 
curve.  As  y  is  decreased,  the  step  here  merges  gradually 
with  the  McCumber-Stewart  curve  without  any  abrupt 
changes,  in  contrast  with  what  happens  in  the  first  lobe. 

The  behavior  of  the  oscillation  along  a  vertical  slice  at 
the  peak  of  the  second  lobe  (-ry  =  3.4)  is  shown  in  Fig.  S. 
Here,  the  situation  is  very  similar  to  that  depicted  in  Fig. 
4  except  at  the  very  top  of  the  step,  where,  with  increasing 
y,  the  / 1  component  just  begins  to  emerge  from  the  back¬ 
ground,  as  seen  in  Fig.  3(c). 

This  tendency  continues  more  markedly  in  Fig.  6, 
which  shows  a  vertical  slice  at  if  —  4.6.  Here,  recalling  the 
situation  depicted  in  Fig.  3,  the  step  is  ZFSl-like  near  the 
top  and  FS2-like  near  the  bottom.  As  in  Fig.  3,  the  tran¬ 
sition  here  between  the  two  oscillation  configurations  is 
relatively  smooth. 

The  information  contained  in  Figs.  2—6  is  summarized 
in  a  highly  schematic  and  approximate  fashion  in  the  in¬ 
sets  of  Fig.  1.  These  show  trajectories  in  the  x-t  plane, 
with  fluxons  and  antifluxons  indicated  by  solid  lines  and 
localized  plasma  waves  indicated  by  dashed  lines.  The  lo¬ 
cations  in  if  of  the  three  insets  in  Fig.  1  indicate  roughly 
the  regions  where  the  corresponding  dynamic  configura¬ 
tions  are  observed  (we  emphasize  once  again,  however, 
that  the  transitions  between  the  various  configurations  are 
gradual). 

A  rather  different  perspective  on  the  nature  of  the  oscil¬ 
lations  is  obtained  by  taking  horizontal  (constant-y )  slices 
through  the  diffraction  pattern.  Two  such  slices,  at 
y  =0.35  and  0.26,  are  shown  in  Figs.  7  and  8.  Figure  7 
shows  the  magnetic  field  dependence  of  the  frequency  f2 
(proportional  to  the  voltage)  in  the  two  slices,  while  Fig.  8 
shows  the  corresponding  behavior  of  the  power  levels,  P\ 
and  Pi,  at  f\  and  f2.  The  salient  facts  that  emerge  from 
these  two  figures  may  be  summarized  as  follows:  (i)  The 
essential,  overall  dependence  of  the  oscillation  frequency 
on  magnetic  field  is  the  inverse  (qualitatively)  of  that  of 
the  diffraction  pattern;  where  the  height  of  the  step  de¬ 
creases  with  field,  the  frequency  increases,  and  vice  versa. 


FIG.  7.  Second-harmonic  frequency  section  through  diffrac¬ 
tion  pattern  at  (a)  y=0,35  and  (b)  y  =0.26.  Arrows  indicate 
hysteresis. 


FIG.  8.  Power  section  through  diffraction  pattern  at  constant 
y.  Powers  measured  at  jc=0  and  normalized  to  4>] 
=  1.0xl0"10.  (a)  First  harmonic  Pt  at  y=0.35.  (b)  First  har¬ 
monic  P |  at  y=0.26.  (c)  Second  harmonic  Pj  at  y  =  0.35.  (d) 
Second  harmonic  P}  at  y  =  0.26.  Arrows  indicate  hysteresis. 


(ii)  In  the  first  lobe  of  the  diffraction  pattern  the  behavior 
is  quite  regular;  the  frequency  increases  monotonically 
with  field,  and  P\  is  the  dominant  component.  We  note 
parenthetically  here  that  this  frequency  behavior  may  be 
different  for  longer  junctions." 19  (iii)  In  the  left  half  of 
the  second  lobe  the  behavior  is  also  regular;  the  frequency 
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FIG.  9.  Partial  section  through  diffraction  pattern  at 
y  ”0.21.  (a)  Tower  levels  of  first  two  Fourier  voltage  com¬ 
ponents  at  x  —  0,  normalized  to  I.Ox  IO-,°.  (b)  Second- 
harmonic  frequency.  Arrows  indicate  hysteresis. 
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here  decreases  monotonically  with  field,  reaching  a 
minimum  in  correspondence  with  the  peak  of  the  second 
lobe,  and  Pi  is  strongly  dominant  (the  / j  component  lies 
within  the  background  noise),  (iv)  In  the  right  half  of  the 
second  lobe  P\  once  again  becomes  dominant.  The  cross¬ 
over  point,  where  P\  a  Pi,  occurs  at  larger  values  of  -q  for 
decreasing  y.  These  observations  suggest  that  it  might  in 
fact  be  more  appropriate  to  refer  here  to  a  third  lobe, 
rather  than  to  a  half  of  the  second  lobe,  (v)  A  strongly 
hysteretic  behavior  of  both  frequency  and  power  is  ob¬ 
served  in  the  right  half  of  the  second  lobe  (i.e.,  third  lobe). 
This  fact  may  explain  why  experimental  measurements 
(see,  e.g.,  Paternd  and  Nordman20  as  well  as  Fig.  12(a) 
below]  often  display  notable  irregularities  in  this  region. 

Figure  9  shows  a  small  section  of  a  horizontal  slice  at 
y  —  0.21  in  the  region  just  under  the  juncture  point  of  the 
first  and  second  lobes.  Since  y=0.26  and  0.35  both  lie 
above  this  juncture  point,  the  curves  in  Figs.  7  and  8  are 
discontinuous  in  this  region.  The  major  conclusions  to  be 
drawn  from  Fig.  9  are  similar  to  those  drawn  above:  (i) 
The  qualitative  shape  of  the  field  dependence  of  the  oscil¬ 
lation  frequency  is  the  inverse  of  that  of  the  diffraction 
pattern,  (ii)  Hysteresis  is  observed  in  17  regions  where  the 
dominant  mode  is  changing. 

IV.  EXPERIMENTAL  RESULTS 

The  numerical  calculations  discussed  above  were  com¬ 
pared  with  measurements  on  overlap  geometry  Nb— Nb- 
oxide— Pb  junctions  having  parameters  comparable  with 
those  used  in  the  calculations.  Although  the  discussion 
below  is  appropriate  for  the  many  junctions  investigated,21 
two  junctions  were  measured  in  detail.  Junction  no.  S10- 
5-1  has  dimensions  479  X  179  /no2,  maximum  zero- voltage 
current  7c0=  1.40  mA,  and  an  estimated  normalized 
length,  L,  of  2.  Junction  no.  65H7  has  dimensions 
467x67  fim1,  7c0=0.53  mA,  and  L  slightly  less  than  2. 
From  independent  measurements  on  similar  junctions  an 
estimate  of  a  and  0  can  be  made.22  The  estimate  is 
reasonably  consistent  with  the  values  a =0.05  and 
0—0.02  used  in  the  calculations,  although  the  experimen¬ 
tal  values  are  probably  somewhat  smaller.  The  parameter 
values  of  both  junctions  are  such  that  the  fundamental 
soliton  frequency  ft  may  be  detected  with  an  2f-band  re¬ 
ceiver  (8—12  GHz).  Any  radiation  at  /2,  however,  is  out¬ 
side  the  frequency  band  of  the  detector  used.  The  mi¬ 
crowave  receiver  had  an  overall  noise  figure  of  about  8 
dB.  By  using  a  spectrum  analyzer  both  the  power  and  the 
frequency  of  microwave  signals  from  the  junction  could 
be  measured.  Generally,  the  received  power  was  25  dB  or 
less  above  the  physical  noise  level  of  the  receiver.  All  data 
discussed  here  were  taken  at  4.2  K. 

To  investigate  the  fluxon  dynamics  the  following  mea¬ 
surements  were  performed:  (i)  I-V  curves  of  the  steps, 
(ii)  The  magnetic  field  dependence  of  the  maximum 
height  of  the  steps,  (iii)  The  magnetic  field  dependence  of 
the  voltage  of  the  steps  with  the  bias  current  as  a  parame¬ 
ter.  (iv)  The  power  and  frequency  of  the  f\  radiation 
emitted  from  the  junction. 

Figure  10  »howi  the  low-voltage  part  of  the  /-  V  curve 
in  zero  magnetic  field  of  junction  no.  S 10-5-1.  Three 
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FIG.  10.  I-V  curve  of  junction  no.  S10-5-1,  showing  three 
zero-field  steps.  The  dashed  curves  indicate  switching. 

ZFS’s  are  seen.  The  high-voltage  part  omitted  in  Fig.  10 
only  shows  the  usual  increase  in  current  at  the  energy  gap 
voltage.  The  dashed  lines  show  switching  at  the  top  of 
the  supercurrent  and  the  three  ZFS’s.  The  7-  V  curve  of 
the  ZFS’s  is  obtained  by  decreasing  the  current  to  obtain  a 
bias  point  just  below  the  foot  of  the  ZFS’s  and  then  in¬ 
creasing  the  current  again.  This  7-  V  curve  is  typical  of 
the  samples  investigated.  Although  not  shown  with  suffi¬ 
cient  voltage  resolution,  the  shape  of  ZFS1  may  be  com¬ 
pared  with  Fig.  2(b). 


H(arb.  units) 


FIG.  1 1.  Junction  no.  S 10-5- 1 .  (a)  Magnetic  field  dependence 
of  the  maximum  zero-voltage  current,  (b)  Magnetic  field  depen¬ 
dence  of  the  height  of  ZFS1/FS2.  (c)  Voltage  tuning  of 
ZFS1/FS2  corresponding  to  bias  level  A  in  (b).  (d)  Magnetic 
field  dependence  of  the  height  of  ZFS2/FS4.  The  dashed  lines 
in  (b)  and  (d)  show  examples  of  parameter  regions  where  /  */ 1 
radiation  was  observed.  The  dotted  lines  show  parameter  re¬ 
gions  where  f  =f  \  radiation  was  not  observed. 
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Figures  11(a)  and  11(b)  show  the  magnetic  field  depen¬ 
dence  of  the  maximum  zero-voltage  current  and  of 
ZFS1/FS2,  respectively.  Note  that  the  lobe  pattern  of 
ZFS1/FS2  is  in  good  qualitative  agreement  with  the  cal¬ 
culated  one  in  Fig.  1.  In  fact,  the  agreement  here  is  at 
least  semiquantitative:  The  ratio  of  the  zero-field  current 
in  Fig.  11(b)  to  that  in  Fig.  11(a)  is  0.63;  the  correspond¬ 
ing  value  from  Fig.  1  is  0.64  <y  <0.65.  The  extrapolated 
zero  of  the  first  lobe  in  Fig.  11(b)  occurs  at  a  field  value 
approximately  equal  to  that  seen  for  the  extrapolated  zero 
of  the  first  lobe  in  Fig.  11(a)  (this  is  a  generally  observed 
experimental  fact).  In  terms  of  the  normalized  field  tj, 
the  extrapolated  zero  of  the  first  lobe  of  the  zero-voltage 
current  occurs  at  rj  =  2  for  very  long  junctions  (the  exact 
value  for  shorter  junctions,  which  is  the  same  for  the 
overlap  and  the  in-line  geometries,  may  be  calculated 
from  the  theory  of  Owen  and  Scalapino23).  Extrapolating 
to  y= 0,  the  first  lobe  in  Fig.  1  also  yields  77  =  2.  The 
dashed  lines  show  examples  of  parameter  regions  where 
/=/ 1  radiation  (at  approximately  9  GHz)  was  observed. 
In  general,  the  /  =fx  radiation  was  observed  in  the  first 
and  third  lobes  but  not  in  the  second  lobe.  This  is  in 
agreement  with  the  results  of  the  numerical  calculations, 
in  particular  Fig.  8.  By  comparison  with  the  calculations 
the  reason  for  the  absence  of  radiation  in  the  second  lobe 
(FS2)  is  that  here  the  /  — /2  radiation  is  at  presumably 
—  18  GHz,  outside  the  range  of  the  receiver.  Figure  1 1(c) 
shows,  for  a  bias  current  corresponding  to  A  in  Fig.  1 1(b), 
the  magnetic  field  tuning  of  the  voltage  of  ZFS1/FS2  (the 
absolute  value  of  the  voltage  is  approximately  35  p\). 
Figure  11(c)  is  in  good  qualitative  agreement  with  the  cal¬ 
culation  of  Fig.  7,  showing  both  the  increase  in  frequency 
(voltage)  as  the  border  regions  of  the  lobes  are  ap¬ 
proached,  and  the  hysteresis  there.  From  the  voltage 
curve,  however,  it  cannot  be  decided  whether  it  is  the 
f  —f  \  or  the  f  —fi  radiation  that  is  dominant.  Thus, 
determination  of  the  fluxon-mode  configuration  requires  a 
measurement  of  the  frequency  of  the  emitted  microwave 
radiation.  Such  a  measurement  is  described  below  (Fig. 
12).  Figure  11(d)  shows  the  lobe  pattern  of  the  second 
step  although  no  corresponding  numerical  calculations 
have  been  performed  for  this  step.  The  dashed  lines  show 
examples  of  parameters  where  f  =f\  radiation  was  ob¬ 
served.  In  general,  no  such  radiation  was  observed  in  the 
first  and  third  lobes;  however,  f  —f  \  radiation  was  ob¬ 
served  everywhere  in  the  second  lobe.  Most  likely,  the 
soliton  configuration  in  the  first  lobe  is  that  of  a  sym¬ 
metric  fluxon-antifluxon  mode  (with  /  =y2 )-  This  is  in 
apparent  contradiction  with  measurements  on  other  sam¬ 
ples4  where  the  f  =f\  radiation  was  also  measured  on 
ZFS2  in  zero  magnetic  field,  but  the  difference  may  sim¬ 
ply  be  a  manifestation  of  the  two  soliton  configurations, 
symmetric  (/=/j)  and  bunched  (/  =/|),  that  have  been 
demonstrated  numerically  for  ZFS2.9  On  junction  no. 
65H7,  in  fact,  the  bunched  mode  was  observed  on  ZFS2. 
In  the  second  lobe  (FS4)  a  three-fluxon— one-antifluxon 
mode  would  give  the  right  frequency  (/()  and  voltage; 
however,  other  configurations  are  also  possible. 

Figure  12(a)  shows  the  diffraction  pattern  for 
ZFS1/FS2  of  junction  no.  65H7.  Qualitatively,  it  is  quite 
similar  to  Fig.  11(b);  however,  the  right  half  of  the  second 
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H  (orb.  units) 

FIG.  12.  Junction  no.  65H7.  (a)  Magnetic  field  dependence 
of  the  height  of  ZFS1/FS2.  The  dashed  lines  show  examples  of 
parameters  where  /  =/i  radiation  was  observed.  The  dotted 
lines  show  parameter  regions  where  /  —f  \  variation  was  not  ob¬ 
served.  (b)  Frequency  of  the  /=/ ,  radiation  corresponding  to 
bias  levels  A,B,C  in  (a). 


lobe  appears  somewhat  anomalous.  For  this  measurement 
the  hysteresis  phenomena  between  lobes  were  observed  but 
not  carefully  mapped.  As  indicated  by  the  lines  A,B,C 
for  different  constant  bias  currents,  f  —f\  radiation  was 
observed  in  the  first  and  third  lobes,  but  not  in  the  second. 
Figure  1 2(b)  shows  a  measurement  of  the  frequency  of  the 
emitted  /=/i  radiation  corresponding  to  bias  currents 
A,B,C  in  Fig.  12(a).  The  positive  frequency  tuning  in  the 
first  lobe  is  in  good  agreement  with  Fig.  7.  The  absence 
of  the  f  —f  \  radiation  in  the  second  lobe  and  the  reap¬ 
pearance  of  such  radiation  at  the  transition  between  the 
second  and  third  lobes  are  also  in  qualitative  agreement 
with  the  numerical  calculation  of  Fig.  8.  As  noted  in  con¬ 
nection  with  Figs.  7  and  9  the  frequency  tuning  resembles 
the  inverse  of  the  diffraction  pattern.  Qualitatively,  curve 
C  in  Fig.  12(b)  is  also  in  agreement  with  the  voltage  tun¬ 
ing  in  the  first  and  third  lobes  of  Fig.  1 1(c). 

V.  CONCLUSIONS 

The  intermediate  regime  L  >  1,  where  the  standard 
analytical  methods  (cavity-mode  theory  for  L  <  I  and  sol¬ 
iton  perturbation  theory  for  L  » 1)  a  priori  do  not  apply, 
was  investigated  numerically  and  experimentally.  A  com¬ 
parison  between  the  experiments  and  the  numerical  calcu¬ 
lations  showed  a  very  good  qualitative  agreement.  Based 
on  this  comparison  it  was  possible  to  identify  the  various 
soliton  modes  in  the  magnetic  field  lobes  of  ZFS1/FS2. 
The  extension  of  single— cavity-mode  theory  to  mul¬ 
timodes  due  to  Enpuku  el  al.1  gave  satisfactory  predic¬ 
tions  of  the  diffraction  pattern  for  ZFS1  and  FS2.  Thus, 
with  caution  elements  from  both  types  of  theories  are  us¬ 
able,  however,  numerical  simulation  is  necessary  for  a 
wider  understanding  of  experimental  observations. 
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calculated  linewidths  with  experimental  results  shows  good  agreement. 
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I.  INTRODUCTION 


Josephson-junction  fluxon  oscillators  continue  to  attract 
research  interest  both  theoretically,  in  studies  of  nonlinear 
wave  dynamics,  and  experimentally,  where  the  very  nar¬ 
row  linewidth  of  the  emitted  microwave  radiation  prom¬ 
ises  potentially  interesting  applications.1  This  very  nar¬ 
row  linewidth  makes  the  numerical  study  of  the  detailed 
dynamics  of  such  oscillators  very  CPU  time  consuming. 
In  order  to  overcome  these  difficulties  we  have  developed 
a  pseudospectral  algorithm  for  solving  the  perturbed 
sine-Gordon  equation  which  describes  the  oscillator.  This 
algorithm  employs  a  Fourier  transformation  of  the  spatial 
variable  together  with  a  finite-difference  approximation  to 
the  time  variable.  The  extensive  use  of  fast  Fourier 
transforms  in  the  algorithm  has  made  the  implementation 
natural  on  a  CRAY-l-S  vector  processor.  The  Fourier 
treatment  of  the  space  variable  requires  spatial  periodicity 
in  the  model.  In  physical  terms  this  means  that  we  are 
studying  a  circular  junction  oscillator  of  the  type  first 
proposed  by  McLaughlin  and  Scott.2  This  device,  as  well 
as  providing  a  convenient  mathematical  model  because  of 
I  periodic  boundary  conditions,  has  in  recent  years  begun  to 
attract  research  interest  in  its  own  right.3-4 

The  paper  is  structured  as  follows.  In  Sec.  II  we 
describe  the  mathematical  model  of  the  circular  junction. 
Details  of  the  numerical  techniques  employed  are  present¬ 
ed  in  Sec.  III.  In  Sec.  VI  we  study  the  behavior  of  the  os¬ 
cillator  under  the  influence  of  a  sinusoidal  driving  term  in 
the  bias  current,  which  models  external  microwave  irradi¬ 
ation.  Section  V  contains  calculations  of  the  linewidth 
under  the  influence  of  Gaussian  white  noise,  which 
models  internal  thermal  noise  in  the  junction.  In  Sec.  VI 
we  compare  our  results  with  existing  experimental  obser¬ 
vations.  In  all  of  the  sections  we  are  focusing  on  a  config¬ 
uration  with  a  single  propagating  fluxon,  which  corre¬ 
sponds  to  the  first  zero-field  step  in  the  current-voltage 
characteristic  of  the  oscillator. 

11.  MATHEMATICAL  MODEL 


As  a  model  for  the  Josephson  tunnel  junction  of  overlap 
geometry  we  use  the  perturbed  sine-Gordon  equation,5 
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1<Pxx=<Ptt—  sinq>  =  a<p,  +  Y +  y(x,t)  .  (2.1) 

•  I  I  i  i 

Here  <p  is  the  quantum  phase  difference  between  the  two 
superconducting  layers  in  the  junction.  Space  and  time 
are  normalized  to  the  Josephson  penetration  length 
\j  —(<t>o/2nj0Ll>)[/1,  and  the  inverse  of  the  plasma  fre¬ 
quency  tap  =  (2trj0/<PuC)w\  respectively,  where  <1>U  is  the 
magnetic  flux  quantum  given  by  <t>0=h/2e  =2.064 
X10-15  Wb.  Lp  and  C  are  the  inductance  and  the  capa¬ 
citance  per  unit  length  of  the  junction.  The  first  of  the 
perturbation  terms  on  the  right-hand  side  of  Eq.  (2.1) 
represents  the  loss  due  to  tunneling  of  normal  electrons,  in 
normalized  units  a  =  G/copC,  where  G_l  is  an  effective 
normal  resistance  per  unit  length.  The  second  term  is  the 
normalized  bias  current  y  measured  in  units  of  j0  the 
maximum  Josephson  current  per  unit  length.  In  this  pa¬ 
per  we  include  a  third  term  7 ;(x,f)  representing  either  an 
externally  applied  sinusoidal  driving  term  connected  to 
the  bias,  or  an  internal  thermal  noise  term  connected  to 
the  loss.  In  this  second  case  we  assume  a  distributed 
Gaussian  white  noise  with  zero  mean  value. 

The  normalized  length  of  the  Josephson  junction 
l~L/\j  is  assumed  to  be  large  compared  with  unity  and 
the  normalized  width  u>=  IV/kj  small  compared  with  un¬ 
ity,  allowing  us  to  use  a  1  +  1  dimensional  model.6  be¬ 
cause  the  aim  of  this  investigation  is  to  isolate  the  influ¬ 
ence  of  the  term  7 )(x,t)  on  the  solution  to  Eq.  (2.1)  we 
avoid  phenomena  connected  with  collision  with  junction 
boundaries  by  considering  a  long  annular  junction.  There¬ 
fore,  we  demand  spatial  periodicity  with  period  /  in  the 
two  physical  quantities,  the  voltage  drop  across  the  junc¬ 
tion: 

4>o Mb 

,2-2' 

j  and  the  current  along  the  junction, 
j  I=~hh<P*  .  (2-3) 

j  i.e.,  boundary  conditions 

q>,l0,t)=<p,ll,t)  ,  (2.4a) 

’  p,(0, »)«?>,(/, f) .  (2.4b)' 
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/  The  fluxon  traveling  wave  solution  to  the  unperturbed 
version  of  Eq.  (2.1)  is  given  by7 

^  =  2sin-l[cn(|',/c )]  ,  (2.5) 

with  £  =  <x-i/r)/[*(l-M2)l/2].  Here  u  is  the  velocity  of 
the  wave  and  k  is  the  modulus  in  the  Jacobian  elliptic 
function.*  Spatial  periodicity  requires  1/(1—  u2),/2 
=  2 nkK(k),  where  n  is  the  winding  number,  i.e.,  the  num¬ 
ber  of  fluxons  minus  the  number  of  antifluxons,  and  K(k ) 
is  the  complete  elliptic  integral  of  the  first  kind.  In  Ref.  9 
it  is  shown  by  Hamiltonian  perturbation  theory2  that  the 
steady-state  fluxon  velocity  dependence  on  the  loss  and 
bias  parameters  is 

u  =  \/(\+(4a'/nY)1)W2  ,  (2.6) 

with  a'  =zaE(k')/k,  where  E(k)  is  the  complete  elliptic 
integral  of  the  second  kind.  For  />  8  (assuming  n  =  1) 
Eq.  (2.5)  reduces  to  the  kink  for  the  infinite  line 
<p  =  4tan-,(e{)  with  £  =  (x  —  ut)/(  1  —  u2)l/2,2  and  the 
velocity  given  by  u  =  1/(1  +(4 a/ny)2]W2.  In  the  numeri¬ 
cal  simulations  we  have  used  /  =  8,  12.8,  20,  and  it  =  1. 
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'  III.  NUMERICAL  TECHNIQUES 

The  very  narrow  Iinewidth  of  the  radiation  emitted 
from  a  Josephson-junction  oscillator  (less  than  1  kHz  at 
10  GHz)10  suggests  that  a  relative  numerical  accuracy  of 
at  least  1 0~7  is  essential.  We  solve  Eq.  (2.1)  numerically 
by  using  a  pscudospcctral  method.11  This  method,  a 
Fourier  transform  treatment  in  space  together  with  a  leap¬ 
frog  scheme  in  time,  has  the  advantage  of  simplicity  and 
high-order  accuracy  in  the  approximations  to  the  space 
derivatives.  Expansion  of  the  fluxon  wave  into  truncated 
series  of  sines  and  cosines  demands  periodicity  not  only  in 
and  tp,  but  also  in  <p  itself.  Observing  that  the  fluxon 
is  a  localized  kink  connecting  two  ground  states  separated 
by  2jt  we  introduce  a  new  periodic  function-  <p  —  2nx/l 
whose  Fourier  representation  we  denote  with  the  su¬ 
perscript  p  =0,  ±  1 . 

Transforming  Eq.  (2.1)  into  the  following  set  of  ordi¬ 
nary  nonlinear  coupled  differential  equations: 

-fc,V(r)— <tf,(t)— F^sin^l 

=a4>?(/)  +  /y8,.o+W'’(f)  ,  (3.1a) 

kp  =  2np/l,  p  =0, ±  1, . . . ,  ±pmlx  (3.1b) 

in  which  Ff  and  Np  are  the  Fourier  components  of  simp 
and  t),  respectively,  and  5p0  denotes  the  Kronecker  sym¬ 
bol,  and  using  second-order  central  differences  to  approxi¬ 
mate  the  time  derivatives  we  get  an  explicit  scheme  for 
the  time  evolution  of  the  Fourier  components 

<t> °j  + ,  =  [  24> °j  -  ( 1  -  a  A  t  /I )  _ , 

-AfVsj’-f  lY  +  Nf)]Al+abt/2) ,  (3.2a) 

4»5+1  =  [(2-At2kp2)<I>J-(l-aA(/2)4>J_, 

-At2(S/+JV/)]/(l-|-aA//2),  |p|>0,  (3.2b) 

\  where  Sf  equals  F?(sin<p)  at  time  jAf,  calculated  each 
Itime  step  by  transforming  (  back  to  x  space,  calculat- 
I  !  I  I  I 
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FIG.  I.  Schematic  diagram  of  numerical  simulation  pro¬ 
cedure. 


ing  simp  and  then  transforming  again  to  k  space  as  indi¬ 
cated  schematically  in  Fig.  1.  ' 

Figure  2  shows  the  computed  <px  as  a  function  of  time 
at  an  arbitrary  point  on  the  junction.  This  signal  consists 
of  an  almost-periodic  sequence  of  pulses.  In  fact,  it  is  the 
deviation  from  perfect  periodicity  that  gives  a  nonzero 
Iinewidth  of  the  radiation.  Since  the  deviation  is  small  it 
is  necessary  to  devise  a  very  accurate  method  for  deter¬ 
mining  the  revolution  periods  T„  for  the  circulating  flux- 
on.  We  do  this  by  calculating  T„  as  the  time  for  the 
mean  value  of  the  phase  over  x  to  change  by  In.  The 
fundamental  frequency  of  the  signal  then  becomes 
/0=l/(r„>,  where  brackets  denote  an  average  value. 
We  take  the  power  spectrum  of  the  signal  near  /0  to  be 
the  distribution  of  the  computer  values  of  1  /T„ . 

Figure  3  shows  the  calculated  Tn's  in  a  computer  ex¬ 
periment  with  the  driving  term  i/  =  0  in  Eq.  (2.1).  As  can 
be  seen  from  Fig.  3,  the  relative  accuracy 
A7Y(r„  >  <  10~8.  In  fact,  examination  of  the  numerical 
output  shows  that  it  is  approximately  7X  10'"<*.  The  long 
transient  arises  from  the  fact  that  the  initial  conditions 
given  by 

<p(x,0)=f(x,0)~ sin-,(y)  ,  (3.3a) 

<p(x,  —  Ar)=/(x,  —  Af)  —  sin-,(y)  ,  (3.3b) 

where  f(x,t)  is  the  fluxon  traveling  wave  solution  to  the 
Unperturbed  sine-Gordon  equation  as  given  by  Eq.  (2.5) 


r  FIG.  2.  Time  dependence  of  the  space  derivative  of  the  flux- 
i:  on  waveform,  showing  the  nth  period  of  revolution  T,  for 
\cr=0.01,  y=0.02,  7=0,  and  /=  12.8. 
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FIG.  3.  Revolution  period  T .  as  a  function  of  revolution 
number  n  for  a=0.01,  >'=0.02,  17=0,  and  /  =  8  showing  high 
level  of  computational  accuracy  achieved, 
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and  sin  (y)  is  the  ground  state,  are  not  exactly  equal  to 
the  final  propagating  configuration. 

We  note  at  this  point  that  the  accuracy  of  the  results 
was  checked  by  doubling  pma,  in  Eq.  (3.1b),  in  order  to 
ensure  that  no  spurious  Fourier  modes  due  to  the  discreti¬ 
zation  in  x  space  are  produced,  and  halving  A /  in  Eqs. 
(3.2).  The  values  used  for  pm„  ranged  from  64  to  256  and 
those  for  A t  from  0.075  to  0.0025,  depending  on  the  pa¬ 
rameters  /  and  y. 

The  computer  program  was  implemented  on  an  IDM 
3033  in  double  precision  (approximately  16  significant 
digits)  and  on  a  CRAY-1  vector  processor  in  single  pre¬ 
cision  (approxximately  15  significant  digits)  using  optim¬ 
izing  FORTRAN  compilers.  In  the  former  case  we  have 
used  the  IMSL-routine  FFT2C  for  fast  Fourier 
transform.12  In  the  latter  case,  by  making  full  use  of  vcc- 
torization  of  the  computer  code  and  the  CRAY  routines 
for  Fourier  transform  and  vector  copying  CFFT2  (Ref. 
13)  and  CCOPY  (Ref.  14)  we  gained  a  speed-up  factor  in 
computing  time  of  22.  Each  long  simulation  requires  typ¬ 
ically  5x  105  time  steps  on  a  512-point  spatial  lattice  and 
uses  approximately  10  min  of  CPU  time  on  the  CRAY- 
1-S  as  opposed  to  approximately  4  h  on  a  scalar  machine. 

Finally,  we  have  compared  the  steady-state  fluxon  ve¬ 
locity,  given  by  u=//(7'II),  with  the  predicted  value 
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from  Hamiltonian  perturbation  theory,  Eq.  (2.6).  The  re¬ 
sult  is  seen  in  Fig.  4.  The  deviation  for  large  bias  values  is 
expected  because  the  perturbation  theory  is  only  valid  for 
small  y  values. 

IV.  SINUSOIDAL  DRIVING  TERM 

In  this  section  we  investigate  the  behavior  of  the  fluxon 
velocity  when  the  driving  term  is  given  by 

-  z>(x,t)  =  T>(t)  =  rj0sin(ft/)  ,  (4.1) 

.  as  a  function  of  the  driving  frequency  ft.  This  might  be 
/  considered  as  a  model  of  microwave  irradiation  of  the 
i  junction.  Using  the  definition  of  the  normalized  momen¬ 
tum 


i  f ' 

pU)=-T  Jo<px<p,dx  , 


and  separating  the  phase  into  a  kink  part  and  a  back¬ 
ground  part15  tp(x,t)=<pk{x,t  )+<p'“(l ),  and  assuming  that 
the  length  of  the  junction  is  large,  allowing  expressions  for 
the  infinite  junction  to  be  used,  we  get  the  following  equa¬ 
tion  for  the  momentum  pk  of  the  kink, 
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FIG.  4.  Difference  between  average  propagation  velocity  as 
computed  numerically  and  calculated  from  perturbation 
theory  u ^  from  Eq.  (2.6)  at  a  function  of  the  bias  for  a=0.0l, 
ij=0,  and  /  =  8. 
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/  FIG.  5.  Revolution  frequency  /,  as  a  function  of  revolution 
number  n  for  sinusoidal  drive,  t>(t)«qgsin(ftr),  with  a-0,01, 
y=0.02,  ft=0.86,  q0=0.0l,  and  /  =  8.  (a)  Numerical  simula¬ 
tion.  (b)  Kink  model. 
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l  x'  dt 

\ 


+ap*=j  y  +  t,0sin(n/)+a^r  +  ^- 


Thus,  the  background  motion  becomes  an  effective 
driving  term  for  the  kink  part.  From  Eq.  (2.1)  we  derive 
the  linearized  equation  for  tp  “=ip“-Fsin-1(y),  assuming 
that  $>  *° « 1, 

+a~^T — +(l-y2)l/J$>“  =  —  i7oSin(ft/)  .  (4.4) 

dt 1  dt 

Combining  Eqs.  (4.3)  and  (4.4)  we  obtain  for  the  kink 
momentum 


pk(t)=  ~  %-  +  — r — ^r-T-^sin(flr—  0X) 

y  4  a  (a2  +  ft2)'/2 


_ Vcfl _ 

|((i-F2)-n2]2+a2n2j 


- - r-r - z — ;  ,  ~  COS(ftt  -02)  , 


Revolution  frequency  f„ 


\/  0,  =  tan-'(n/a) 


02  =  tan-||aft/[(l-y2),/2-ft2]|  . 

The  instantaneous  kink  velocity  is  then  calculated  from 
pk  —  u/(  1  _MJ)I/2.  In  order  to  compare  this  approximate 
theoretical  description  with  the  numerical  result  we  calcu¬ 
late  the  nth  period  Tn  according  to  the  formula 

r  *•  - 1 +  ,  . 

I.  ud t=l , 


'  •  u  =p*/[l  -f  (p*)2]1/2  . 

i 

I ■  Figures  5—7  show  a  comparison  of  the  results  from  this 

linearized  model  and  from  numerical  simulations  of  Eq. 
(2.1)  with  ft  =0.86,  0.89,  and  1.10,  respectively.  In  all 
cases  it  is  seen  that  the  kink  model  is  able  to  reproduce 
the  fluctuations  in  the  revolution  frequency  f„  =  \/TK  in 
.great  detail. 

\  As  a  measure  of  the  amplitude  of  the  frequency  flue- 

r-f — h - ; -  ' 

'  i  Revolution  frequency  fn 
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I  FIG.  6.  Revolution  frequency  /.  as  a  function  of  revolution 
number  n  for  sinusoidal  drive,  »j( / )  =  »joSin( fit ),  with  a=0.01, 
ly-0.02,  n-0.89,  tjo«*O.OI,  and  /* 8.  (a)  Numerical  simula¬ 
tion.  (b)  Kink  model. 
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I  FIG.  7.  Revolution  frequency  f„  as  a  function  of  revolution 
I  number  n  for  sinusoidal  drive,  Tf<r)~wsin (ft/),  with  a»0.0l, 
y«0.02,  ftol.io,  T]am0.0\,  and  /■> 8.  (a)  Numerical  simula¬ 
tion.  (b)  Kink  model. 
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/Hamiltonian  perturbation  theory  for  the  fluctuations  Au 
in  the  fluxon  velocity  leads  to  the  power  spectrum  for 
::  au,10 


jo2n-i4)’/2 


<u2  +  a1 
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I  .  '  Driving  frequency  ft 

FIG.  8.  Standard  deviation  of  revolution  frequency  a/  as  a 
function  of  driving  frequency  ft.  Solid  curve,  numerical  simula¬ 
tion;  dashed  curve,  kink  model;  parameters,  a=0.01,  y  =0.02, 
tj0=0.01,  and  7  =  8. 
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tuation,  which  is  essentially  the  linewidth  of  the  oscillator, 
we  have  calculated  the  standard  deviation  of  the  revolu¬ 
tion  frequency  o  f  —  —  </«))2)],/2  for  values  of  the 

cyclic  driving  frequency  ft  between  0.4  and  2.0. 

The  full  curve  in  Fig.  8  shows  the  results  from  the  nu¬ 
merical  simulation  and  the  dashed  curve  those  from  the 
kink  model.  The  kink  model  predicts  a  resonance  just 
below  the  plasma  frequency  ft=  1,  whereas  the  numerical 
simulation  yields  this  peak  at  a  somewhat  lower  frequen¬ 
cy.  Moreover,  the  numerical  results  exhibit  a  hysteresis 
not  seen  in  those  of  the  kink  model  and  a  difference  in 
scale.  The  discrepancy  in  resonance  frequency  and  hys¬ 
teresis  behavior  is  attributable  to  the  fact  that  we  have 
used  a  linearized  kink  model.  Presumably,  the  use  of  a 
higher-order  expansion  in  Eq.  (4.4)  would  yield  a  behavior 
analogous  to  that  of  a  soft  nonlinear  spring16  thus  reduc¬ 
ing  these  discrepancies.  It  is  not  clear,  however,  to  what 
extent  the  difference  in  scale  would  be  resolved  by  such  a 
refinement. 

V.  GAUSSIAN  WHITE  NOISE 


with  the  average  velocity  u0  given  by  Eq.  (2.6).  By  a 
'  Fourier  transform  of  Eq.  (5.4)  we  obtain  the  autocorrela- 
tion  function  for  Au  as  an  exponential 

„  ,  .  a«i‘  l  uo> 


*  Thus  Au(f)  is  a  normal  process  with  zero  mean  and 
standard  deviation 11 


oJl-uJ)5 


OAu  =  ‘ 


'  Defining  the  period  of  a  fluxon  revolution  according  to 
Eq.  (4.6)  we  calculate  the  average  frequency  fluctuation  as 
:the  average  of  the  instantaneous  frequency  fluctuation 
'Au  //  over  one  average  period  of  revolution 

;  '  Y~~rt+(r) 

A/=7~y /,  Au//dr  .  (5.7) 

From  Eq.  (5.7)  it  follows  that  A/  has  a  normal  distribu- 
'  tion  with  zero  mean  and  the  standard  deviation,19 

1  2u0  1  —  exp(  —  al/u0)  1/2 

^A/=yOAy  —  1  U0 - — -  .  (5.8) 


The  term  ij(x,r)  in  Eq.  (2.1)  is  here  considered  to  be  \ 
Gaussian  white  noise  with  zero  mean  (ij(x,f))=0  and  !:■; j 
-  autocorrelation  function  :-.]4  • 

:;i;: 

Jtn{g,r)=(Tr(x,/)ij(x+g,f  +  r))=afa£)5lr) .  (5.1)  y..; 

—  The  variance  of  the  noise  a J  is  connected  with  the  loss  a  y 

and  the  absolute  temperature  T  through17  jfe 

_  a\-AirakT ,  (5.2) 

|V  t 

where  k  is  the  Boltzmann  constant.  .  j 

In  the  vectorized  algorithm  we  find  it  convenient  to  in- 

—  troduce  the  noise  term  in  p-t  space,  Nf{t),  as  f 


Distribution  density 
3.0x10s  1 


2.0x10s 


1.0x10s 


7V'(r)=F-'(<71,exp(i(0,+0J)|  ,  (5.3)  j. 

where  F-1  denotes  the  Fourier  transform  from  a  to  t  s' 
space,  and  Qf  and  0W  are  stochastic  variables  uniformly  j.. 
distributed  between  0  and  2ir,  with  an  upper  limit  in  p  ’ 
and  a  of  pmil »  1/2 Ax  and  cumal»rr/At,  Ax  and  A t  being 
the  resolution  in  space  and  time,  respectively.  Standard 

'l  !  i  I  !  i 
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A  numerical  simulation  with  o,  =  8.8x  10~4  is  seen  in 
Fig.  9  showing  a  typical  frequency  distribution  of  A / 
about  the  fundamental  frequency  /0  =  u0//.  The  connec- 
;•  tion  between  the  standard  deviation  and  the  half-power 
linewidth  is 

A/,/2  =  v/8h^2aA/  (5.9) 

j  when  A/  is  normal  distributed. 

\  Figures  10  and  11  show  a  comparison  of  the  standard 
/deviation  predicted  by  this  model  Eq.  (5.8)  and  the  results 
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|  FIO.  9.  Distribution  of  revolution  frequency  /,.  Numerical 
!  simulation  with  Gaussian  noise  drive:  a—0.01,  y— 0.034, 
!<7,«8.8x10-\  and  /-20. 
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FIG.  10.  Standard  deviation  of  revolution  frequency  oy  for 
white  Gaussian  noise  drive  as  a  function  of  bias  current  y,  for 
a  =  0.0I,  /  =  8,  and  tr,  =  0.01,  0.05,  0.10,  and  0.20.  (a)  Numeri¬ 
cal  simulation,  (b)  Hamiltonian  perturbation  theory. 
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from  the  numerical  simulations  for  the  lengths  /= 8  and 
/  =  20,  respectively.  As  can  be  seen,  the  model  is  able  to 
predict  the  right  qualitative  dependence  on  the  length,  the 
noise  amplitude,  and  the  bias,  but  the  model  predicts  an 
overall  standard  deviation  that  is  about  a  factor  of  10  too 
large.  The  reason  for  this  discrepancy  is  at  present  not 
known. 

In  closing,  we  note  that  for  y  values  near  0.3  it  was 
found  necessary  to  augment  the  time  resolution  (by  reduc¬ 
ing  At)  to  avoid  spurious  peaks  in  Fig.  10(a).  The  ex¬ 
istence  of  such  spurious  peaks  might  be  an  indicator  of 
the  onset  of  chaotic  behavior  at  nearby  points  in  parame¬ 
ter  space.  In  fact,  parameter  values  y  =  0.3  with  a=0.01 
lead  to  chaotic  creation  of  fluxon-antifluxon  pairs  in  the 
study  reported  by  Eilbeck.20 

VI.  COMPARISON  WITH  EXPERIMENTS 

The  rapidly  decreasing  linewidth  with  increasing  bias 
shown  in  Figs.  10  and  11  is  in  qualitative  agreement  with 
the  experimental  observations  of  Fig.  1  in  Ref.  10. 

To  compare  quantitatively  the  calculated  results  with 
these  experiments  we  use  in  Eq.  (5.2)  data  reported  by 
i  Scott  et  al .5  For  the  junction  No.  N25L,  assuming  a  tem- 
\  perature  of  4  K,  Eq.  (5.2)  gives  a„  =0.0052.  Noting  from 
Fig.  9(a)  that  Of  scales  linearly  with  an,  we  calculate  from 
Eq.  (5.9)  a  normalized  half-power  linewidth 

A/|/2  =  5.5X  10-7  at  y=0.2.  Taking  as  the  normalized 
resonance  frequency  /0  =  w0//~ 0. 125  we  calculate  a  rela¬ 
tive  linewidth  A/|/2//0=4.4X  10~6.  The  physical  reso¬ 
nance  frequency  for  junction  No.  N25L  was  2.3  GHz.5 
This  yields  a  physical  linewidth  of  10  kHz.  Comparing 
:  with  the  experimental  results  shown  in  Fig.  1  of  Ref.  10 
:  and  noting  that  y=0.2  corresponds  to  a  bias  point  near 
f  the  bottom  of  the  zero-field  step,  we  find  excellent  agree- 
:  ment.  The  same  calculations  for  junction  No.  N53C,5 
again  for  T= 4  K  and  y=0.2,  yield  A/|,j//0 
=  2.3X  10-6.  The  physical  resonance  frequency  for  junc¬ 
tion  No.  N53C  was  8.3  GHz,  whicn  leads  to  a  physical 
:  linewidth  of  18  kHz,  once  again  in  excellent  agreement 
\with  experimental  results. 


VII.  CONCLUSIONS 

Computational  studies  of  the  linewidth  of  the  radiation 
emitted  by  Josephson  junctions  require  extremely  high 
resolution.  For  this  reason  we  developed  a  pseudospectral 
method  for  solving  the  nonlinear  dynamical  equation 
describing  a  circular  Josephson  junction  oscillator.  Be¬ 
cause  the  algorithm  makes  heavy  use  of  fast  Fourier 
transforms  it  was  implemented  on  a  CRAY-1  vector  pro¬ 
cessor.  Driving  terms  corresponding  to  physically  realis¬ 
tic  situations,  i.e.,  sinusoidal  microwave  irradiation  and 
internal  thermal  noise,  were  considered.  In  the  second 
case  the  computational  results  were  compared  with  exper¬ 
imental  results  reported  in  the  literature,  and  excellent 
qualitative  and  quantitative  agreement  was  found.  In  ad¬ 
dition,  in  both  cases  we  have  compared  the  computational 
results  with  approximate  analytic  results  based  on  pertur¬ 
bation  theory.  Here  the  agreement  was  qualitatively  good, 
but  quantitative  discrepancies  were  found,  indicating  a 
need  for  further  development  of  perturbation  theory. 
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Superconducting  quantum  interference  devices  (SQUIDs)  respond  chaotically  to  external  oscillating  fluxes.  The  small 
deviations  from  one-dimensional  return  maps  are  investigated.  Within  a  narrow  region  of  parameter  space  a  sequence  of  period 
doublings,  windows  with  odd  periods  and  chaotic  behaviour,  intermitlency  and  bifurcation  between  coexisting  attractors  of  low 
dimension  are  found. 


A  periodically  driven  rf  superconducting  quantum 
interference  device  (SQUID)  consisting  of  a  ring  with 
a  single  Josephson  junction  has  been  investigated  re¬ 
cently  both  by  means  of  analogue  circuits  [1]  and 
computationally  (2 ,3] .  Like  many  other  physical  sys¬ 
tems  the  device  exhibits  deterministic  chaos  (intrinsic 
noise).  In  this  letter  we  investigate  in  detail  a  particu¬ 
lar  portion  of  parameter  space  which  exhibits  ex¬ 
tremely  rich  details  of  the  dynamics.  Furthermore, 
we  demonstrate  that  the  choice  of  attractor  (in  cases 
with  more  than  one  attractor)  may  depend  on  the  ini¬ 
tial  conditions,  i.e.  there  can  be  coexisting  attractors. 

The  order  parameter  for  the  Josephson  junction,  <t>, 
satisfies  the  differential  equation  [2] 

<p"  +  e<p'  +  sin  0  =  a(y  sin  coD  t  -  4>)  ,  (1) 

where  primes  signify  differentiation  with  respect  to 
the  dimensionless  time  t  given  by  t  =  7XC0o/2jryc)1/2. 
Here  T is  laboratory  time,  <>0  =  hj 2\e\  is  the  flux 
quantum,  C  and  Jc  are  the  capacitance  and  the  maxi¬ 
mal  critical  current,  respectively.  The  loss  parameter, 
e,  is  given  by  e2  =  0o/2nC/c/?2  where  R  is  the  resis¬ 
tance  of  the  weak  link.  The  external  flux  is  assumed 
to  be  sinusoidal  with  frequency  and  amplitude  y. 
In  the  additional  term,  a<p,  a  =  <Pq/2tiJcL,  where  L  is 
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the  inductance  of  the  ring.  Due  to  the  presence  of 
this  term,  which  corresponds  to  a  confining  quadratic 
potential,  less  chaos  might  be  expected  in  this  system 
than  in  the  single  Josephson  junction  without  the  a<t> 
term  j2).  We  have  shown  that  eq.  (1)  does  not  possess 
the  Painlev6  property  (4)  in  accordance  with  the  fact 
that  chaos  in  the  system  does  occur  for  certain  initial 
conditions.  The  proof  is  obtained  by  converting  eq. 
(1)  into  four  non-linear  coupled  autonomous  first- 
order  differential  equations  through  differentiations 
and  trigonometric  substitutions  [5].  The  resulting 
system  can  be  shown  to  possess  other  movable  singu¬ 
larities  than  poles.  Therefore  the  system,  and  thus  eq. 
( 1),  is  not  of  Painlevg  type. 

Eq.  (1)  was  integrated  numerically  by  means  of 
IMSL  routine  [6]  DVERK-I.  Our  results  are  shown  as 
the  return  maps  in  figs.  1  and  2.  As  in  ref.  (2]  we  plot 
X{n  +  1)  versus  X(n)  =  +  A)  with  tn  given  as  the 

nth  zero  crossing  of  <p'  and  the  constant  A  arbitrarily 
chosen  as  A  =  0.134.  The  numerical  simulations 
showed  that  the  transients  had  died  out  before  n  = 
300.  Iterations  from  n  =  300  to  n  =  600  are  included. 
In  view  of  the  fixed  rf  driving  our  general  return  map 
construction  could  be  substituted  by,  e.g.,  a  simple 
strobing  at  the  driving  period.  Differing  constructions 
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should  not  affect  the  topological  information  carried 
by  the  map.  The  computational  procedure  proposed 
by  Henon  [7]  was  used.  Of  the  four  parameters  in  cq. 

( 1 )  e,  i  and  cjd  were  kept  constant  throughout  at  the 
values  chosen  in  fig.  3  of  ref.  (2| :  e  =  0.5,  7=10,  and 
u>D  =  0.5. 

In  fig.  1  the  parameter  a  was  varied  from  a  = 
0.16435  to  a  =  0.173.  In  the  first  run,  a  =  0.16435, 
initial  conditions  0(0)  =  <fi'( 0)  =  0  were  used  and  a 
period-2  solution  was  found  (as  indicated  in  fig.  2  of 
ref.  [2[ ).  In  the  following  run,  a  =  0.1644,  the  final 
values  of  0  and  <p'  from  the  previous  run  were  used  as 
initial  conditions.  This  procedure  was  used  through¬ 
out  in  a  scries  of  computer  experiments  with  increas¬ 
ing  values  of  the  parameter  a.  However,  for  each  value 
of  a,  eq.  (1)  was  also  solved  for  initial  conditions  0(0) 
=  0'(O)  =  O. 

When  the  resulting  return  maps  (for  n  =  300-600) 
differed  the  map  is  included  in  fig.  2.  For  a  =  0.1644 
a  return  map  indicating  low-dimensional  chaotic  be¬ 
haviour  was  found  (for  both  sets  of  initial  conditions). 
As  the  parameter  a  was  raised  to  a  =  0. 1 65  only  small 
changes  in  the  return  map  occurred.  However,  at  a  = 
0.167  we  observed  a  shift  from  chaotic  behaviour  into 
a  period-7  solution.  At  a  =  0.169  the  return  map  again 
exhibits  chaotic  behaviour,  this  time  as  a  period-3 
window.  At  o  =  0.1695  the  periodic  behaviour  is  re¬ 
established  as  a  period-3  solution  which,  however, 
vanishes  again  at  a  =  0.17  where  we  find  the  same  dou¬ 
ble  humped  chaotic  return  map  as  in  fig.  3c  of  ref. 

[2] .  Note,  however,  that  if  initial  conditions  0(0)  = 
0'(O)  =  0  are  used  chaos  already  occurs  for  a  =  0.1695 
as  shown  in  fig.  2.  Already  for  a  =  0. 1 7007  the  dou¬ 
ble  humped  return  map  depopulates  into  a  chaotic 
period-5  window,  which  has  already  vanished  again 
for  a  =  0.1705  (fig.  1).  The  double  humped  envelope 
curve  for  this  value  is  very  similar  to  the  curve  ob¬ 
tained  for  a  =  0.17.  A  closer  analysis  of  the  return 
map  shows  that  tire  system  stays  for  a  number  of 
iterations  on  the  lower  hump  (which  has  a  double¬ 
curve  structure).  Then  it  moves  to  the  upper  hump 
(single-curve  structure)  for  a  number  of  iterations  be¬ 
fore  it  goes  back  to  the  lower  hump.  This  intermitten- 
cy  cycle  continues  as  far  as  we  have  followed  it. 

When  the  parameter  a  is  raised  to  the  value  a  = 
0.1706  two  things  may  happen:  The  system  may 
choose  the  attractor  which  corresponds  to  the  lower 
hump  of  the  curve  (as  shown  in  fig.  1)  or  the  attractor 
which  corresponds  to  the  upper  hump  of  the  curve 


(as  shown  in  fig.  2).  Thus  a  special  type  of  bifurcation 
occurs  between  a  =  0. 1705  and  0. 1706.  By  means  of 
computer  experiments  we  have  checked  that  the 
choice  depends  on  whether  the  initial  conditions  for 
the  run  already  lie  on  the  lower  attractor  or  the  upper 
attractor,  respectively.  Also,  we  have  found  that  ini¬ 
tial  conditions  0(0)  =  0'(O)  =  0  lead  always  to  the  up¬ 
per  attractor  as  shown  in  fig.  2.  The  following  runs 
demonstrate  that  the  system  stays  on  the  preferred 
attractor  as  or  is  increased.  For  u  =  0.17125  the  lower 
attractor  (fig.  1)  and  the  upper  attractor  (fig.  2)  have 
depopulated  into  period-2  windows.  For  or  =  0.17129 
period- 16  solutions  (on  the  attractors)  occur  in  both 
cases.  For  or  =  0.17  1 3,  0. 17  14,  and  0. 17 15  a  sequence 
of  period-8,  period-4,  and  pcriod-2  solutions  is  ob¬ 
served  similarly.  For  a  =  0.173  the  period- 1  solution 
indicated  in  fig.  2  of  ref.  |2)  is  found. 

The  chaotic  behaviour  indicated  by  “c”  at  e  =  0.5 
and  a  =  0.16-0.17  in  fig.  2  of  ref.  [2  J  has  an  extreme¬ 
ly  detailed  structure  within  a  small  portion  of  param¬ 
eter  space.  In  particular,  note  that  the  “single  hump” 
return  map  reported  in  ref.  [2]  is  in  fact  extremely 
structured.  Since  the  phase  space  is  two-dimensional  a 
purely  one-dimensional  return  map  is  not  expected, 
the  closeness  depending  sensitively  on  the  parameter 
e.  Similar  small  deviations  from  a  one-dimensional 
circle  map  have  been  found  in  ref.  [8J  in  a  dc  +  ac 
driven  single  pendulum. 
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